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INTRODUCTION 


Any  active  space  system  will  generate  a  self-induced  environment  due  to  the  release  of  various 
effluents  and  the  effluents'  interacdcHis  with  the  ambient  space  environment.  The  pemiibed  local 
environment  may  interfere  negatively  with  space  system  operations,  and  it  is  therefore  important 
to  know  what  environmental  perturbations  to  expect,  and  predict  the  effect  of  these  perturbations 
cm  space  system  operaticms.  Space  systems  include  space  platforms  such  as  Space  Station 
FreedcHn  (SSF)  and  Strategic  Defense  Initiative  (SDI)  type  weapons'  platforms,  satellites,  orbital 
transfer  vehicles  (OTVs),  the  Shuttle  Orbiter,  and  future  space  transportation  systems.  The 
altitude  regime  we  have  focused  cm  is  low  Earth  orbit  (LEO),  roughly  spanning  the  altitudes  from 
200  to  1000  km. 

The  main  effluent  source  is  the  operaticm  of  various  thrusters,  which  include  attimde  control  and 
reboost  operations,  as  well  as  venting  from  extrmely  high  mass  flow  rate  chemical  power 
systems.  Even  though  the  thruster  plume  exhaust  is  initially  neutral,  it  can  become  ionized 
through  a  series  of  interacticxis  with  the  ambient  envimunent.  Such  icmization  mechanisms 
include  chaise  exchange  with  the  ambient  oxygen  ions,  i^ioto-ionization  due  to  the  solar  extreme 
ultra  violet  (EUV)  photcxi  flux,  and  Critical  Icxiization  Velocity  (CTV)  processes.  So  in  addition 
to  an  ambient  background  plasma,  there  may  be  a  local  plasma  enhancement  due  to  ionization  of 
the  neutral  effluents.  The  perturbed  plasma  provides  an  artificial  envimunent  for  any  high- 
voluge  surface  (ui  the  spacecraft,  such  as  solar  arrays.  The  high-voltage  surfaces  may  be  exposed 
to  the  space  environment  either  through  design  or  by  erosion  of  insulation.  For  certain  plasma 
and  neutral  density  regimes,  arcing  and  electric  breakdown  may  be  triggered  at  the  high-voltage 
surfaces,  causing  sputtering  and  erositxi  of  surfaces  such  as  protective  and  thermal  coatings.  The 
sputtered  material  can  redeposit  on  other  sensitive  surfaces  such  as  solar  arrays,  radiators, 
sensors,  and  windows.  Erosion  of  thennal  coatings  will  affect  the  global  thermal  management  of 
the  space  system,  and  the  loss  of  protective  coatinp  exposes  the  surfaces  to  both  harmful  ultra 
violet  (UV)  radiation  from  the  Sun  and  to  the  chemically  active  ambient  atomic  oxygen,  which 
can  wear  down  surfaces  by  reacting  with  them.  All  these  processes  reduce  the  efficiency  and 
lifetimes  of  space  systems,  and  must  be  understood  in  order  to  eflectively  design  a  complex 
space  system. 

In  this  effort  we  have  concentrated  on  modeling  the  self-induced  environments  expected  at  high- 
voltage  solar  arrays.  This  macroscopic  environmeru  serves  as  input  to  a  microscopic  model  of  die 
arcing  process,  developed  by  Professor  Daniel  Hastings  and  Dr.  Mengu  Cho  at  the  Massachusetts 
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Institute  of  Technology  (MIT).  The  goal  is  to  provide  a  predictive  capability  for  the  occurrence 
and  prevention  of  arcing  on  high-voltage  surfaces. 

The  contract  objective  expressed  in  the  Statement  of  Work  is  to "  develop  a  model  of  the  neutral 
and  charged  particle  self-induced  envirtmment  around  a  space  platform  /  vehicle".  To  keep  the 
analysis  relevant,  we  have  focused  on  two  realistic  effluent  scenarios  for  space  systems  in  LEO: 
(1)  Massive  H2  release  ( kg/s  mass  flow  rate)  from  an  SDI-type  chemical  power  system  (e.  g.. 
Neutral  Particle  Beam  Platform),  and  (2)  SSF  hydrazine  thrusters. 

After  the  structure  geometry  and  thruster  coniiguratiai  were  defined,  the  neutral  effluent 
flowfields  were  modeled.  The  effluent  density  regimes  range  fnxn  continuum  inside  and 
immediately  outside  the  thrusters,  via  transitional  to  rarefied  flow  beyond  that  The  continuum 
part  of  the  flowfields  were  modeled  with  a  Method  of  Qiaracteristics  (MOC)  code  under 
ccxnpany  Internal  Research  and  Develt^ent  (IRAD).  The  rarefied  pan  of  the  plumes  were 
characterized  with  direct  simulatitxi  Monte  Carlo  (DSMQ  methods,  developed  by  Dr.  G.  A. 

Bird.  The  ionization  of  the  neutral  plumes  was  determined  by  solving  a  set  of  continuity 
equations  for  the  neutral  and  plasma  species. 

The  total  calculated  neutral  and  plasma  environment  was  then  coupled  with  the  microscopic 
arcing  model  developed  at  MIT  to  predict  arcing  probability  and  characteristics  for  the  modeled 
space  system. 

RESULTS 

Neutral  Piaritmlion 

Two  realistic  effluent  releases  have  been  modeled;  one  is  the  venting  out  of  a  high  mass  flow  rate 
SDI-type  chemical  power  system,  and  the  other  is  the  operation  of  varicms  hydrazine  thrusters  in 
the  SSF  environment 

Chemical  Power  System  Exhaust 

One  of  die  concepts  within  the  SDI  program  is  a  space  based  neutral  particle  beam  weapons 
platform  (see  figure  l).Such  a  device  would  »;celerate  negative  hydrogen  ions  through  an 
electron  stripper,  hence  generating  a  beam  of  neutral  hydrogen.  Power  levels  on  the  order  of 
megawatts  are  required  to  accelerate  the  hydrogen,  and  one  concept  is  powered  by  a  chemical 
power  system  fueled  by  hydrogen  and  oxygen.  The  combusdtm  products  are  water  and  hydrogen, 
where  die  water  is  contained  in  an  internal  cycle,  but  the  hydrogen  is  released  to  the  ambient 
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icHiosphere.  The  hydrogen  gas  would  be  vented  through  a  series  of  nozzles  at  mass  flow  rates  on 
the  order  kg/s  for  several  minutes.  The  hydrogen  molecules  will  become  partly  ionized  after  exit 
and  perturb  the  local  plasma  environment.  The  perturbed  environment  can  interact  with  the  beam 
itself  (e.  g.,  attenuadtm,  deflection),  and  with  high-voltage  surfaces  on  the  space  platform. 

We  have  modeled  such  a  release  through  a  superscmic  nozzle  at  a  mass  flow  rate  of  0.83  kg/s. 

The  continuum  part  of  the  flow  was  modeled  by  an  MOC  code,  and  the  transition  line  between 
continuum  and  rarefied  flow  was  defined  at  a  Bird  number  of  0.04.  The  Bird  number,  which  is  a 
condition  for  transition  fitmi  continuum  to  rarefied  flow,  is  discussed  in  more  detail  in  the 
DSMC  overview  in  Appendix  A. 

A  schematic  view  of  the  nozzle  and  spacecraft  is  showed  in  figure  2.  The  nozzle  radius  is  28  cm 
and  the  nozzle  half  angle  is  25  degrees.  For  the  case  we  present  here,  the  spacecraft  radius  was  50 
cm,  and  the  high-voltage  solar  array  surfaces  would  be  located  around  the  structure  upstream  of 
the  nozzle  exit  plane.  The  amdnuum  prt^rties  alcmg  the  Bird  transition  line  were  input  to  a  2-D 
DSMC  code.  Ibe  DSMC  omcept  is  described  in  more  detail  in  Appendix  A.  Required  inputs  are 
composititm,  molecular  number  density,  velocities,  flow  angles,  and  temperature.  The  DSMC 
rurm  were  executed  on  a  Staixlent  super  mini-computer,  and  run  times  were  on  the  ord^r  of  15-20 
CPU  hours. 

The  H2  ctnistant  number  density  contours  ate  shown  in  figure  3.  For  this  case,  the  thruster  was 
assumed  to  be  venting  into  a  vacuum.  The  DSMC  code  does  have  provisions  for  including  the 
effects  of  an  ambient  free  stream.  The  flee  stream  has  little  effect  on  high  density  thruster  plumes 
(like  here),  and  to  keep  the  analysis  general  and  independent  of  altitude  and  ambient  additions,  a 
vacuum  boundary  omdititxi  was  implemented.  The  nozzle  is  located  to  the  right  altxig  the  z-axis, 
and  tlw  number  densities  vary  ffran  lO^i  m*^  at  the  nozzle  lip  to  10^®  m'^  at  the  arrays 
upstream  of  the  nozzle  exit  Typical  ambient  neutral  densities  in  LEO  are  of  the  order  10^^  to 
10^^  m*^.  Figure  4  depicts  the  static  pressure  atrxmd  the  space  system,  and  these  range  from 
10~^  to  1(H  torr  at  the  structure  walls.  This  is  greater  than  the  ambient  pressure  of  t^roximately 
lO"®  to  l(f^  torr. 

SSF  Hydrazine  Thrusters 

The  SSF  leboost  and  attitude  control  will  be  provided  by  25  Ibf  (-1 1 1  N)  and  55  Ibf  (~245  N) 
monopropellant  hydrazine  (N2H4)  thrusters.  In  «idition,  the  SSF  envimunent  will  be  affected 
by  the  Shuttle's  25  Ibf  bipn^llant  hydrazine  Vernier  Reactitxi  Control  System  (VRCS)  thrusters 
during  Shuttle  docking  and  berthing  maneuvers.  The  bipn^Uants  are  fueled  by  MMH-N2O4, 
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where  MMH  stands  for  monomethyl  hydrazine  (CH5N2)  and  N2O4  is  the  oxidizer.  We  have 
simulated  both  mcmoprop  and  biprop  activity  for  a  series  of  operating  and  geometrical  parameters 
including  various  mass  flow  rates  and  locations  with  respect  to  high  voltage  surfaces.  The  SSF 
will  be  assembled  over  a  total  of  17  Shuttle  tnissicms.  and  the  station  orbital  attitude  will  vary.  In 
the  final  configuration  (figure  S).  the  SSF  will  orbit  with  the  arrays  facing  the  ram  and  the  truss 
perpendicular  to  the  orbital  velocity  vector.  During  the  first  assembly  phases  shown  in  figures  6 
and  7,  however,  the  station  will  be  orbiting  in  a  so-called  arrow  mode,  with  the  truss  oriented 
along  the  velocity  vector.  To  maintain  orbital  altitude,  the  station  will  need  to  reboost  a  few 
times  per  year  during  the  arrow  mode  phase.  The  reboost  will  be  provided  by  25  Ibf  monoprops 
located  tm  propulsitxi  modules,  which  will  be  operating  either  over  and  across  the  solar  arrays,  or 
almig  the  truss  in  the  tqrposite  direction.  Figure  8  shows  a  scenario  with  the  thrusters  operating 
across  one  bank  of  solar  arrays,  with  the  approxim^ly  10  by  30  m  arrays  located  about  S  m 
downstream  of  the  thruster  exit  plane.  The  reboost  (^rations,  which  will  last  a  few  hours  every 
few  months,  will  generate  a  strongly  perturbed  local  neutral  /  plasma  environment  at  the  high- 
voltage  (160  V)  solar  arrays.  In  our  analysis  we  have  computed  the  expected  envirorunent  to 
determine  if  arcing  on  the  arrays  is  a  concern. 


The  thrusters  that  were  simulated  have  an  exit  radius  of  1.62S  iit  (4.13  cm)  and  an  expansion 
ratio  of  100:1.  To  demcxistrate  the  characteristics  of  biprop  versus  monoprop  hydrazine  thruster 
flowfields,  we  have  also  simulated  SSF  25  Ibf  hydrazine  thrusters  (grating  on  bipropellant 
MMH-N2O4,  even  though  this  fuel  mix  is  not  used  on  the  SSF  (only  on  the  Shuttle).  The 
theoretical  exhaust  composidtnis  for  the  mmoprcq)  arxl  biprop  hydrazine  thrusters  are  shown  in 
the  following  table  1: 


Table  1 :  Theoretical  Exhaust  Composition  of  Hydrazine  Thrusters 

Monopropeliant  Hydrazine  (N2H4) 

Exhaust  Specie 

Molar  Fraction 

Hydrogen  (H2) 

0.48 

Nitrogen  (N2) 

0.29 

Ammonia  (NH3) 

0.22 

Bipropellant  Hydn 

izine  (MMH-N204) 

Exhaust  Specie 

Molar  Fraction 

Hydrogen  (H2) 

0.20 

Nitrogen  (N2) 

0.31 

Water  (H20) 

0.32 

Carbon  Monoxide  (CO) 

0.11 

Carbon  Dioxide  (C02) 

0.06 
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The  hydrazine  thrusters  are  located  on  propulsion  modules  that  are  replaced  when  the  iuel  is 
expended.  Towards  the  end  of  the  module  lifetime,  the  mass  flow  rate  decreases,  and  to  simulate 
this  development,  thrust  levels  of  17, 9,  and  2.S  Ibf  were  modeled  in  addition  to  the  nominal  25 
Ibf.  Only  the  25  and  2.5  Ibf  flovrfields  will  be  presented  here;  the  17  and  9  Ibf  levels  were 
presented  in  the  annual  technical  repon  of  31  July,  1991. 

The  thruster  locations  with  respect  to  the  solar  arrays  are  illustrated  in  figure  8.  The  arrays  will 
constantly  be  tracking  the  Sun,  however,  during  reboost,  the  arrays  will  be  feathered  to  some 
degree  to  avoid  having  the  thrusters  firing  directly  at  the  arrays.  The  arrays  rotate  both  in  alpha 
and  beta  directions,  and  the  local  neutral  densities  at  the  surface  of  the  arrays  will  be  a  function 
of  array  attitude.  The  case  that  was  looked  at  here  has  the  array  located  in  a  horizontal  position 
with  an  angle  alpha  of  90  degrees,  and  an  angle  beta  of  0  degrees. 

The  SSF  is  traveling  at  an  orbital  velocity  of  ^roximately  7.7  km/s.  This  is  expressed  in  the 
DSMC  simulation  as  a  free  stream  of  atomic  oxygen  flowing  past  the  arrays  along  the  y-axis  in 
figure  8.  At  LEO  altitudes,  auxnic  oxygen  is  the  predominant  neutral  species  with  a  number 
density  of  about  10*  - 10^  cm*3.  The  nominal  temperature  of  the  arrays  varies  between  199  and 
333  K,  depending  cm  the  solar  heating;  most  of  the  simulations  were  done  for  a  nominal  array 
temperature  of  250  K.  For  most  of  the  cases  full  timrmal  acccmunodaticm  at  the  array  surface  has 
been  assumed.  The  effect  of  varying  the  acccmunodaticm  coefficient  is  addressed  later  on. 

25  Ibf  Monopropeliant  Hydrazine  Thruster 

This  is  the  tiiruster  that  will  actually  be  used  to  provide  reboost  for  the  SSF  during  the  arrow 
mcxle.  The  ccmtinuum  part  of  the  flowfield  was  calculated  using  MOC  methcxls  starting  inside 
the  thruster.  The  flow  properties  (compositicm,  density,  vector  velocity,  temperature)  along  a  line 
separating  continuum  and  transiticmal  flow  were  used  as  input  to  the  DSMC  simulation.  Figure  9 
shows  the  ccmstant  number  density  ccmtours  for  this  thruster.  The  ambient  free  stream  is  fnxn  left 
to  right,  so  the  thruster  is  operating  into  the  wake  of  the  SSF.  The  simulated  area  is  axi- 
symmetric  with  the  thruster  center  line  being  the  symmetry  axis.  The  ccmtinuum  part  of  the  flow 
iirunediately  outside  the  thruster  exit  shows  up  as  a  white  regicm,  and  the  white  lines  and  streaks 
in  the  color  plots  only  signify  various  computational  regions  in  the  DSMC  simulaticm;  the 
se|  araticm  between  regicms  is  ncx  physical,  but  purely  a  result  of  the  way  the  color  grai^cs 
software  interpolates  the  data  (or  rather  the  lack  of  interpolaticm).  In  the  simulaticm  itself,  the 
regions  are  ccmtinuous. 


5 


The  exhaust  plume  initially  exhibits  a  nominal  expwsion  and  density  drop,  but  downstream, 
close  to  the  solar  anay  (which  is  oriented  out  of  the  paper  plane),  a  density  "pile-up"  is  evident, 
causing  the  local  neutral  molecular  number  density  to  reach  levels  of  lO^^  m'^,  which  is  six 
orders  of  magnitude  greater  than  the  ambient  density.  Note  that  this  is  the  total  neutral  density. 
The  monoprop  hydrazine  exhaust  consists  of  hydrogen,  nitrogen,  and  ammonia,  and  the 
relatively  lighter  hydrogen  will  show  a  greater  expansion  into  the  backflow  region  than  the  other 
two  species.  The  DSMC  code  does  keep  track  of  the  various  species,  but  for  our  analysis,  we 
have  just  assumed  the  species  molar  fraction  to  be  fixed  at  all  regions  of  the  flowfield.  Figures  10 
and  1 1  show  the  static  pressure  and  the  temperature  for  the  exhaust  plume.  The  static  pressure  is 
just  p  =  1/2  pv2,  where  p  is  the  mass  density,  and  V  is  the  velocity  of  the  fiow.  The  static 
pressure  can  be  seen  to  first  decrease  out  of  the  tK>zzle  to  levels  of  10*^  torr,  before  building  up 
towards  the  array,  reaching  values  of  0.01  torr.  This  compares  to  ambient  pressures  of  about  10‘^ 
torr.  The  temperature  is  really  an  average  temperamre  averaged  over  the  translational,  rotational, 
and  vibrational  temperatures  of  tl^e  molecules.  For  this  particular  case,  full  thermal 
accommodation  was  assumed,  which  means  that  the  flow  incident  to  the  array  surface  is  brought 
to  the  array  surface  temperamre.  The  flow  initially  cools  off  as  it  expands  out  of  the  thruster 
nozzle,  but  as  the  flow  slows  down  and  stagnates  towards  the  array  surface,  it  heats  up,  due  to  the 
conversion  of  molecular  translational  kinetic  energy  to  internal  enei^gy.  The  flow  heats  up  to 
about  17CX)  K  in  a  layer  along  the  array,  before  it  cools  off  again,  this  time  being  accommodated 
to  the  array  surface  temperamre  fixed  at  250  K.  The  net  heat  flux  to  the  solar  array  surface  varies 
from  26  W/m^  (3  m  downstream  of  the  thruster)  to  191  W/m^  (8  m  downstream  of  the  thruster). 

Figure  12  shows  the  absolute  exhaust  velocity  for  an  array  temperamre  of  333  K.  TTie  flow 
initially  expands  out  of  the  nozzle  at  a  speed  of  approximately  3  km/s.  The  orbital  velocity  of  the 
SSF  is  7.7  km/s,  which  means  that  the  exhaust  travels  at  7.7  -  3  =  4.7  km/s  with  respect  to  the 
ambient  free  stream.  The  flow  slows  down  as  it  approaches  the  array,  becoming  stagnant  (with 
re^a  to  the  array)  about  2  m  downstream  of  the  nozzle  exit  plane.  Figures  13  and  14  show  the 
axial  and  radial  components  of  the  flow,  respectively.  Right  outside  the  thruster  ex-^  plane,  4  m 
off  the  symmetry  axis,  one  can  see  how  the  axial  velocity  component  is  reversed,  the 

exhaust  to  expand  into  the  backflow  at  speeds  as  high  as  0.8  km/s;  this  will  impact  the  CIV 
iotuzatitm  in  this  region,  as  will  be  seen  later.  The  radial  velocity  component  is  zero  along  the 
thruster  axis,  and  about  3  km/s  off  the  nozzle  lip.  Towards  the  array  the  axial  cennponent  is 
suppressed  again,  showing  reversal  and  negative  axial  velocity  close  to  the  surface. 
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25  Ibf  Bipropellant  MMH>N204  Thruster 

This  next  case  is  a  simulation  of  an  SSF  thruster  fueled  by  bipropellant  hydrazine,  which  is  what 
the  Shuttle  VRCS  thrusters  run  on.  The  flowlleld  number  densities  shown  in  figure  IS  resemble 
the  ones  for  the  monoprop,  showing  a  density  "pUe-up”  at  the  array  8  m  downstream  of  the 
thruster  exit  plane.  The  colored  streaks  in  this  flowfield  are  attributed  to  numerical  non-physical 
phenomena.  The  exhaust  number  densities  for  the  biprop  are  generally  lower  than  for  the 
monoprop  by  almost  half  an  order  of  magrutude.  This  can  be  attributed  to  the  difference  in  the 
mass  of  the  exhaust  species;  the  biprop  mixture  has  a  lower  mass  per  mole  than  the  monoprop, 
which  causes  it  to  expand  faster  and  which  ultimately  leads  to  a  less  dense  plume. 

2,5  Ibf  Monopropellant  Hydrazine  Thruster 

This  is  basically  a  SSF  thruster  ruiming  at  a  mass  flow  rate  one  order  of  magnitude  less  than 
nominal.  Figure  16  shows  the  constant  number  density  contours  for  this  case;  the  solar  array 
inhibits  the  free  expansion  of  the  thruster  plume,  but  there  is  no  substantial  density  build-up  at 
the  array  as  demonstrated  for  the  25  Ibf  case.  The  flow  temperatures  are  displayed  in  figure  17 
and  they  are,  as  one  would  except,  significantly  lower  than  for  the  25  Ibf  case.  A  one  order  of 
magnitude  decrease  in  the  mass  flow  rate  causes  the  temperamre  to  be  cut  in  half,  reaching 
approximately  750  K  in  a  layer  7  m  downstream  of  the  thruster  exit  plane. 

2,5  Ibf  Bipropellant  MMH*N204  Thruster 

Hgure  18  shows  the  constant  number  density  contours  for  this  case.  Again,  the  overall  flow 
density  is  reduced  compared  to  the  25  Ibf  case.  By  cranparing  with  the  2.5  Ibf  monoprop 
densities  shown  in  figure  16,  we  see  the  effea  of  the  lower  plume  molecular  mass,  causing  the 
lighter  biprop  plume  to  expand  faster  and  become  more  rarefied. 

Surface  Accommodation 

The  effect  of  varying  the  array  surface  accommodaticm  was  also  looked  at  The  concept  of 
surface  accommodation  is  still  not  completely  understood,  and  any  conservative  analysis  should 
check  the  effect  the  accommodation  coefficient  on  the  flow  properties. 

A  gas  in  contaa  with  a  solid  surface  is  one  of  the  most  frequently  occurring  boundary  conditions. 
Unfortunately,  it  is  also  one  of  the  least  well  understood.  Despite  a  great  deal  of  theoretical  and 
experimental  study,  the  simple  models  of  specular  and  diffuse  reflecticMi  that  were  originally  put 
forward  by  Maxwell  in  1879  remain  the  most  generally  useful  models  for  practical  applications 
[Ref.  1]. 
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Specular  reflection  is  perfectly  elastic  with  the  molecular  velocity  component  normal  to  the 
surface  being  reversed,  while  that  parallel  to  the  surface  remains  unchanged.  This  is  a  useful 
reference  state  for  analytical  studies,  but  does  not  occur  for  real  gas-solid  combinations.  With  this 
model,  there  is  no  thennal  or  viscous  boundary  layer  adjacent  to  a  smooth  solid  surface. 

In  diffuse  reflection  the  velocity  of  each  molecule  after  reflection  is  independent  of  its  incident 
velocity.  However,  the  velocities  of  the  reflected  molecules  as  a  whole  are  distributed  according 
to  the  half-range  equilibrium  Maxwellian  distribution  for  the  molecules  that  are  directed  away 
from  the  surface.  This  distribution  is  for  the  temperamre  T,.  which  may  differ  from  the 
temperature  Ty^f  of  the  surface.  The  extent  to  which  the  reflected  molecules  have  their  temperature 
adjusted  toward  that  of  the  surface  is  indicated  by  the  accommodation  coefflcient  ,  defined  by 


(li-Qw 

Here,  and  q,.  are  respectively  the  incident  and  reflected  energy  fluxes,  while  17^  is  the  energy 
flux  that  would  be  carried  away  in  diffuse  reflecticMi  with  The  range  of  is  from  zero 

for  no  accommodation  to  unity  for  complete  thennal  acconunodation. 

Experiments  with  engineering  surfaces  in  contact  with  gases  at  temperatures  on  the  order  of  the 
standard  temperature  indicate  that  the  reflection  process  approximates  diffuse  reflection  with 
complete  thennal  accommodation.  This  may  be  a  consequence  of  such  surfaces  being 
microscopically  rough  with  the  incident  molecule  suffering  multiple  scattering,  or  of  the 
molecules  being  mcRnentarily  trapped  or  adsorbed  on  the  surface  (e.  g.,  if  the  molecules  react 
chemically  vnth  the  surface).  On  the  other  hand,  experiments  with  carefully  prepared  and  cleaned 
surfaces  indicate  that  the  accommodatitm  coefficient  can  be  significantly  less  than  one.  As  the 
translational  energy  of  the  incident  molecules  relative  to  the  surface  becomes  large  in  comparison 
with  the  value  corresponding  to  the  surface  temperature,  complete  thermal  accommodaticm 
becomes  less  readily  achieved. 

In  figure  1 1  we  saw  the  flow  temperature  for  full  thermal  accommodation,  i.  e.,  an 
accommodation  coefficient  of  one,  at  the  array.  The  solar  array  surface  is  at  250  degrees  K, 
whereas  the  flow  temperature  is  much  hotter,  reaching  peaks  of  approximately  1700  K.  The 
molecules  hitting  the  surface  are  cooled  down  to  the  surface  temperature  (250  K)  -  this  is  evident 
in  figure  1 1  where  there  is  a  relatively  cool  temperature  layer  next  to  the  array.  For  a  specular 
reflection  case,  the  molecules  hit  and  leave  the  surface  without  changing  temperature,  and  the  gas 
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heats  up  towards  the  surface  as  it  is  being  slowed  down,  converting  kinetic  energy  to  internal 
enei^y.  This  is  demonstrated  in  figure  19. 

The  effect  of  thennal  accmnmodadtxi  cn  the  local  density  is  displayed  in  figures  9  and  20. 
According  to  the  perfect  gas  law,  pV  =  nRT,  when  the  temperamre  drops,  the  density  increases  - 
and  vice  versa.  So  for  the  case  with  complete  thermal  accommodation,  the  relatively  lower 
temperature  along  the  array  leads  to  a  higher  local  density.  The  density  "pile-up"  in  figure  9  can 
thus  be  seen  to  be  about  a  faaor  of  two  greater  than  in  figure  20.  Subsequendy,  if  the  local 
density  at  the  array  was  barely  enough  to  trigger  arcing  cm  the  high-voltage  array,  a  factor  of  two 
decrease  in  density  could  help  keep  the  array  from  arcing.  In  that  sense,  it  would  be  preferable  to 
have  a  surface  of  thermal  accommodation  coefficient  zero,  which  is  a  surface  exhibiting  specular 
reflection.  It  should  also  be  noted  that  for  full  thermal  accommodation,  varying  the  array 
temperature  from  199  to  333  K  (min  and  max  due  to  solar  heating),  had  no  detectable  effea  on 
the  local  flow  number  density. 

Ways  to  keep  a  surface  as  specularly  reflecting  as  possible  are: 

(1)  Maintain  smooth  surface. 

•  One  could  presumably  use  (Mily  surfaces  that  are  smooth,  however,  one  must  also 
remember  that  in  space  surfaces  may  be  subject  to  micrometeoroid  impacts,  chemical 
erosion,  and  UV  deformation  •  which  all  tend  to  make  the  surface  pitted  and  non-uniform. 

(2)  Chemically  inactive  or  inert  surfaces. 

•  This  obviously  depends  on  the  incident  molecules  too.  Atomic  oxygen  is  very  reactive  and 
would  be  hard  to  keep  inactive.  However,  the  solar  array  is  made  up  of  solar  cells  with 
cover  glass  and  protective  coating,  metallic  interconnects,  and  supporting  structure.  Both 
cojpptT  and  aluminum  interconnects  rapidly  form  oxides  with  the  attnnic  oxygen,  and  this 
should  halt  further  chemical  activity  involving  oxygen  cm  the  surface.  The  cover  glass 
coatings  also  form  oxides. 

However,  for  the  case  with  the  thruster  operating  altmg  the  array,  the  incident  mcdecules  ate 
ptedomirumtly  the  exhaust  species,  which  are  H2,  N2,  and  NH3  for  monoprcq)  hydrazine,  and 
H2.  N2,  H2O,  CO,  and  CO2  for  biprop  MMH-N2O4.  Neither  N2  rror  H2  are  very  reactive,  which 
should  indicate  a  greater  degree  of  specular  reflection  than  for  a  pure  attnnic  oxygen 
environment  Water  and  ammonia  have  chemical  reactivities  somewhere  between  atomic  oxygen 
and  nitTOgen  /  hydrogot. 
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In  conclusi(»i.  if  the  spacecraft  designer  expects  density  conditions  that  are  borderline  arc 
triggering,  it  is  worthwhile  to  employ  speculary  reflecting  materials.  Research  and  experiments  to 
determine  the  acctmunodatitm  coefficients  for  relevant  materials  would  be  needed  in  order  to 
pick  optimum  materials.  Thermal  acammodation  furthermore  causes  drag,  which  is  another 
reason  for  studying  the  issue.  I  am  not  aware  of  any  such  optimization  taking  place  in  the 
community  today. 

Ionization  of  the  Neutral  Effluents 

Having  computed  the  neutral  flowfield  distributitHis.  we  can  now  turn  on  the  various  possible 
ionization  processes.  lonizadCHi  of  the  neutral  effluents  will  perturb  the  local  plasma  distribution 
and  possibly  aflea  arcing  thresholds  and  arcing  frequencies  at  the  high-voltage  solar  arrays. 

The  neutrals  can  becone  itxiized  through  a  series  of  interactions  with  the  ambient  environment. 
The  most  important  ionization  mechanisms  are: 

Ion-Molecule  Reactions 

These  reactions  lead  to  the  exchange  of  charge  between  species.  Charge  exchange  reactions 
include  charge  transfer  processes  in  which  me  or  more  electrms  are  transferred,  as  well  as  im- 
atrni  or  im-molecule  interchange  (charged  rearrangement)  processes  in  which  a  charged  or 
neutral  atomic  or  molecular  species  is  transferred  similar  to  an  ordinary  chemical  reactim  [Ref. 

2].  Molecules  such  as  H2,  H2O,  and  CO2  react  with  ionosi^ric  oxygen  ions  in  a  charge  transfer 
process  as  follows  [Ref.  3]: 


XY  +  0+  XO+  +  Y 

where  XY  is  the  released  molecule,  and  XO'*'  is  a  positive  molecular  im  which  can  rapidly 
recombine  with  electrons.  Dissociative  recmtbinatim  is  represented  by  the  reaction 


XO+  +  e*  X*  +  O* 


where  X*  and  O*  are  electronically  excited  states  which  will  radiate  visible  light  There  is  of 
course  also  straightforward  recombinatim  such  as: 

XO+  +  e*  XO 

An  example  of  an  im-molecule  interchange  is  [Ref.  4]: 
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OH+  +  H2  H2O+  +  H. 


There  is  no  momentum  exchange  or  sharing  of  kinetic  energy  in  the  charge  transfer  process. 
Furthermore,  charge  exchange  does  not  provide  net  ionization  gain  in  a  global  sense,  however,  if 
a  fast  (7.7  km/s  orbital  velocity)  ambient  plasma  is  replaced  by  a  slower  (effluent  exit 
velocity  ~  3  km/s)  local  plasma  of  ionized  effluents,  the  local  plasma  density  will  increase. 

Photo-Ionization 

A  solar  UV  phoum  of  energy  ht),  where  h  is  the  Planck  constant  and  v  is  the  photon  frequency, 
can  be  absorbed  by  an  atom  or  molecule,  exciting  it  to  a  higher  state,  or  photo-ionizing  it  when 
ho  is  greater  than  the  ionization  energy  [Ref.  S]: 

N2  +  hu  (<  796  Angstrom)  N2'*'  +  e‘ 

Dissociative  photo-ionizatitHi  is  an  additional  source  of  atomic  ions  [Ref.  6]; 

N2  +  ho(<  510  A)  ->  N+  +  N  +  e* 

In  die  latter  case,  the  phottms  must  have  enough  energy  to  both  ionize  and  dissociate  the 
molecule.  The  ion  production  rate,  q,  depends  cm  the  number  density  of  the  ionizable 
constituent,  n,  its  icmization  cross  section  for  the  particular  photon  wavelength  interval,  o,  and 
the  local  photon  flux, 


qson^ 

The  solar  EUV  flux  as  a  fvmction  of  waveloigth  is  well  documented.  The  icmizing  phoum  flux 
decreases  from  its  value  outside  tte  atmosphere  as  the  result  of  absorption.  The  qitical  depth, 
T(X),  is  a  parameter  that  specifies  die  attenuation  of  solar  irradiance  by  the  atmo^diere,  and  the 
ion  production  rate,  adjusted  for  absorption,  can  be  expressed  as: 

q  =  on(z)^exp(-x) 

where  z  depicts  the  altitude.  Hence,  the  ionization  rate  decreases  by  a  factor  e  for  an  qiticai  depth 
of  1.  The  ion  production  rate  has  a  maximum  where  dqfdz  =  0,  which  conespcmds  to  the 
condition  x  =  1  [Ref.  2] 
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The  photo-ionization  process  is  usually  slower  than  the  charge  exchange  processes,  and  all 
species  with  a  photo  lifetime  of  10^  s  or  greater  will  in  general  be  depleted  more  rapidly  by 
chemical  or  transport  processes  than  by  photolysis  [Ref.  4],  Species  like  H2  and  CO2  have 
photo-ionization  rates  of  about  10*^  s'^.  which  corresponds  to  photo  lifetimes  of  1/10"^  =  10^  s. 
Shorter  lifetimes  are  observed  for  H2O  (10^  s)  and  Ni(CO)4  (100  s). 

Critical  Ionization  Velocity  (CFV) 

The  concept  of  CIV  was  first  introduced  by  Affvin  [Ref.  7]  in  1954  as  a  component  of  his  theory 
for  the  formation  of  the  solar  system.  As  a  neutral  gas  moves  through  a  magnetized  plasma,  a 
rapid  itmization  of  the  neutrals  takes  place  if  the  kinetic  energy  of  the  neutrals  relative  to  the 
plasma  exceeds  the  ionization  potential  of  the  neutrals.  This  then  defines  a  critical  velocity, 

Vciv> 


^MVciy2  =  e(t) 

where  M  is  the  neutral  mass,  above  which  CTV  ionizaticm  will  take  place.  Although  equating  the 
neutral  kinetic  energy  to  the  ionization  potential  in  the  calculation  of  a  threshold  for  ionization  of 
the  neutrals  may  appear  at  first  glance  to  be  an  obvious  step,  it  is  important  to  keep  in  mind  that 
the  kinetic  energy  needs  to  be  at  least  twice  that  given  in  the  equation  above  for  an  ionizing  ion- 
neutral  collision  to  occur  (as  is  easily  seen  in  the  center  of  mass  frame).  This  is  because  in 
collisions  between  ions  and  neutrals,  only  a  fractitMi  (one-half  for  equal  mass  ion  and  neutral)  of 
the  kinetic  energy  of  the  collisitm  is  available  for  icmization  of  the  neutral  [Ref.  8].  However, 
strong  evidence  has  been  put  forward  in  lab  experiments  that  CIV  ionization  does  indeed  occur 
and  that  the  threshold  is  quite  near  or  only  slightly  above  that  predicted  by  the  above  equation. 

The  QV  process  depends  on  several  steps  [Ref.  9].  Scnne  mechanism,  such  as  charge  exchange 
between  the  fast  neutrals  and  cold  plasma  itms,  or  friroto-ionization  of  a  fraction  of  the  neutrals, 
creates  energetic  ions.  The  relative  motion  between  these  itms  and  the  plasma  electrons  creates  an 
unstable  situation  which  gives  rise  to  the  growth  of  plasma  waves.  These  waves  serve  to  remove 
energy  from  die  fast  ions  and  transfer  it  to  the  electrons.  The  electrtms  in  the  tail  of  the  energy 
distribution  with  energy  higher  than  e^  can  ionize  neutrals.  This  creates  additional  ions  moving 
at  the  neutral  velocity.  These  newborn  ions  feed  energy  into  the  unstable  waves,  which  in  turn 
further  heat  the  elections,  thereby  allowing  the  cyclic  energy  transfer  to  be  sustained.  The 
exponential  increase  in  plasma  density  is  a  cmisequence  of  this  energy  transfer  cycle. 
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Attempts  to  reproduce  the  CIV  effect  in  space  experiments  have  provided  mixed  success,  but 
research  is  still  going  strong  trying  to  demortstrate  the  CIV  process  in  space  and  unravel  the 
process.  The  most  ccMivincing  evidence  that  the  CIV  effect  does  occur  in  a  space  plasma  is  from 
the  Projea  Porcupine  flight  reported  by  Haerendel  [Ref.  10].  A  neutral  barium  (V^y  =  2.7  km/s) 
release  demonstrated  initial  icmization  of  30  ±  10%  of  the  total  neutral  population  that  had  a 
transverse  velocity  greater  than  the  critical  velocity.  Other  ionization  mechanisms  than  CIV  were 
ruled  out  as  insufficient. 


The  velocities  required  to  trigger  CIV  ionization  are  often  too  high  for  the  effect  to  be  important 
in  LEO.  This  is  because  the  critical  velocities  often  exceed  the  orbital  velocity  (7-8  km/s)  of  the 
space  system  and  its  effluents.  However,  by  thrusting  into  the  ram  direction,  the  exhaust  plume 
velocities  with  respect  to  the  ambient  magnetized  plasma  becomes  the  sum  of  the  orbital  and  the 
thruster  exit  velocities  which  may  well  exceed  the  critical  velocity.  Table  2  shows  the  critical 
ionization  velocities  for  a  few  common  exhaust  species: 


Table  2:  Critical  Velocities  for  a  few  Common  Hydrazine  Exhaust  Soecies 

Soecies 

Vhi/  r  km/sl 

H2 

38.5  [Ref.  11] 

N2 

10.3  [Ref.  12] 

NH3 

10.7  [Ref.  12] 

H20 

11.6  [Ref.  12] 

CO 

9.7  [Ref.  11] 

C02 

7.7  [Ref.  121 

For  the  numerical  simulations  of  the  ionization  CIV,  we  have  taken  advantage  of  imtization  rates 
calculated  by  Dr.  Rodger  Biasca,  formerly  of  MIT,  but  currently  at  the  Phillips  Lab/Hanscom  Air 
Force  Base  [Ref.  12,13].  These  rates  were  obtained  from  a  one-dimensional  implicit  Particle-in- 
Cell  (PIQ  simulation  assuming  a  constant  density,  infinity  background  of  neutrals.  The  implicit 
PIC  code  is  diftferent  from  previous  PIC  simulations  reported  in  the  literature  in  that  the  electrtm- 
ion  mass  ratios  are  real.  The  simulations  assume  a  neutral  beam  (corresponding  to  the  thruster 
plume)  propagating  across  a  magnetic  field  in  a  background  plasma.  It  was  found  that  the 
electrai  number  density  increases  as  roughly 

Ng(t)/Ng(0)  =  c''ci^* 

where  is  the  anomalous  or  CIV  ionization  rate.  This  can  be  rewritten  as: 
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Currently,  (Hily  Xe,  Ne,  CO2,  NO,  and  N2  have  been  simulated,  of  which  CO2  and  N2  are 
^plicable  to  the  hydrazine  thruster  plumes  we  are  modeling.  The  ionization  rates  computed  by 
Biasca  have  been  plotted  as  a  function  of  both  neutral  velocity  and  neutral  density,  and  we  have 
determined  numerical  rates  based  cm  functional  fitting  to  the  plots.  Obviously,  no  CIV  icmization 
occurs  at  velocities  below  the  critical  ionization  velocity,  and  the  rates  increase  quadratically  or 
faster  with  increasing  neutral  velocity.  The  r^es  are  also  proportional  to  the  neutral  density  up  to 
a  pdnt  where  the  neutral  number  density  is  about  2  x  10^  times  the  initial  electron  density.  For 
greater  neutral  densities,  the  electron  collisicm  rates  become  on  the  order  of  the  lower  hybrid 
frequency,  in  which  case  it  is  difficult  for  the  electrms  to  develop  the  hi^  energy  tail  on  the 
electrcm  distribution  ftincticMi  that  plays  such  a  big  role  in  CIV.  The  functionally  fitted  CIV  rates 
for  CO2  and  N2  are  of  the  form 

v/n  =  (K  +  n„/nj)(V„/V^^)5 


where  Q  is  the  ion  gyro  frequency,  K  is  a  constant,  and  n^  are  the  neutral  and  electron 
densities,  respectively,  and  V,^  and  V^^y  ^  ^  neutral  velocity  and  the  critical  ionization 
velocity,  respectively.  The  fitted  CIV  rates  are  shown  in  table  3: 


Table  3:  Functionally  Fitted  CIV  Ionization  Rates  (Valid  for  Vn  >  Vriu  only) 

X  =  Vn  /  Vpjv  and  y  =  Nn  /  Nemi 

C02: 

Neutral  Density  Regime 

■Ociv !  Qc02 

1  e6  <  Nn/Nefol  <  2.5e7 

[0.7  +  0.45(y/1  .e6  - 1 )]  (x/1 .5)^ 

2,5e7  <  Nn/Nomi  <  7e7 

[1.4  +  15.6(1  -v/7e7)l  (x/1 .5)5 

N2: 

Neutral  Density  Regime 

Vciv  /  iiN2 

1  e6  <  Np/Nefo)  <  1 .5e7 

[0.33  +  0.5(y/1e6  - 1 )]  (x/1 .5)i> 

1 .5e7  <  Nn/No/m  <  6e7 

10.6  +  7.2(1  -v/6e7)l  (x/1 .5)5 

Note  that  these  rates  are  valid  only  when  the  neutral  velocity  is  equal  to  or  greater  than  the  QV 
velocity!  The  rates  can  be  seen  to  be  proportional  to  the  rreutral  density,  and  a  very  strong 
function  of  the  ratio  between  tseutral  and  CTV  velocity.  Numerical  rates  are  unfortunately  not 
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available  for  CO  and  NH3,  which  with  critical  ionization  velocities  of  9.7  and  10.7  km/s, 
respectively,  could  well  exhibit  CIV  effects  at  orbital  velocities. 

Numerical  Simulation  of  the  Ionization  Process 

In  order  to  simulate  the  itmization  of  the  neutral  thruster  effluents,  we  established  a  set  of  2-D 
continuity  ^uations  for  all  neutral  and  ion  species  in  the  flow.  The  set  of  hyperbolic  equations 
were  of  the  form: 

dn/dt  +  df/dx  +  df^dy  =  S  -  L 

where  n  denotes  the  neutral  or  irni  density,  f  denotes  the  2-D  flux,  and  S  and  L  are  source  and 
loss  terms,  respectively.  This  set  of  ion  and  neutral  cmitinuity  equations  were  solved  by  time 
advancing  the  flux  through  a  computatimal  grid,  using  the  donor-cell  method.  The  solution  is 
first  order  in  both  time  and  space.  In  order  for  the  numerical  solution  to  be  stable,  the  time  step 
must  obey  the  Courant  stability  ccmdititm,  which  is  expressed  as 

|V|At/Ax<l 

where  V  is  the  flux  velocity,  At  is  the  time  step,  and  Ax  is  the  spatial  step. 

The  initial  neutral  densities  and  vector  velocities  were  provided  by  the  DSMC  runs  in  a  semi- 
raiKlom  dau  grid  format.  This  data  needed  to  be  organized  before  it  was  input  to  the  itMuzation 
model.  A  set  of  DI-3000  subroutines,  available  on  our  VAX  work  station,  were  combined  to 
create  an  interpolated  organized  neutral  grid.  The  gridding  algorithm  used  a  triangulation 
scheme,  followed  by  bivariate  interpolation  with  a  fifth-degree  polynomial. 

The  fact  that  we  are  looking  at  plume  induced  effects  occurring  along  the  solar  arrays  means  that 
we  are  interested  in  ionization  occurring  only  on  time  scales  on  the  order  of  the  plume  exhaust 
propagation  time  passed  the  stto  arrays.  Fm*  thruster  exit  velocities  on  the  order  3  km/s,  and  a 
stdar  array  dimension  of  10.S  m,  it  takes  the  rreutral  effluents  about  3.S  milliseconds  to  blow  past 
the  array.  Ionization  htqrpening  after  that  will  not  influence  the  local  environment  at  the  array. 

Some  of  the  features  of  our  ionization  model  are  as  follows: 


I 
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New  icms  created  as  a  result  of  ion-molecule  reactions,  photo-ionizadon,  and  CIV  effects 
were  assiuned  to  keep  traveling  at  their  initial  neutral  velocities,  i.  e.,  no  collisional 
momentum  was  exchanged. 

Magnetic  effects  were  omitted,  which  is  acceptable  at  the  dimensions  we  are  looking  at.  The 
itm  gyro  radii  range  frnn  about  20  to  100  m,  and  the  gyro  periods  are  between  20  and  60  ms. 

The  effect  of  the  solar  anay  voltage  (Hi  the  ions  is  not  modeled.  However,  for  negative  array 
voltages  and  positive  i(His,  the  newly  formed  ions  will  be  attracted  by  the  array,  which  will 
cause  the  local  i(Mi  density  to  increase. 

The  ambient  plasma  is  assumed  to  be  O'**.  The  directional  velocity  of  the  O'*'  can  be  varied  to 
simulate  ram  and  wake  effects,  i.  e.,  the  thruster  curating  in  different  directions.  The 
ambient  O'**  density  is  initially  set  to  a  ccHistant  number  everywhere  in  the  numerical  grid; 
once  the  itm-molecule  reacti(His  start,  the  O'**  starts  to  deplete,  but  new  ambient  O'**  is 
supplied  with  the  ambient  free  stream. 

The  electron  density  was  initially  set  e<]ual  to  the  UHal  i(Hi  denaty  after  each  time  step.  This 
implies  diat  there  are  etKMigh  electrons  available  at  all  times  to  neutralize  the  newly  created 
l(x»l  plasma.  This  is  usuaUy  not  a  bad  assumption,  however,  for  cases  when  the  CIV  effect 
really  takes  off  and  the  i(Hiizati(Hi  rates  grow  exponentiaUy  with  time,  (me  must  verify  tha 
the  ambient  electron  flux  is  sufficient  to  keep  the  newly  created  local  plasma  overall  charge 
neutral  a  all  times. 


25  Ibf  Monopropellant  Hydrazine  Thruster 

As  shown  in  table  1.  the  exhaust  plume  for  this  thruster  consists  of  H2,  N2,  and  NH3.  The 
various  ionization  processes  these  species  are  subjea  to,  are  shown  in  table  4: 


Table  4:  Ionization  Rates  Monopropellant  Hydrazine 

Ion-Molecule  Reactions 

Reaction 

Reaction  Rate  fcm3/sl 

0+  +  H2  0H+  +  H 

1.7e-9 

[Ref.  4] 

0H+  +  ©■  ^ 

0*  +  H 

7.5e-8  (300/Te)’/2 

[Ref.  4] 

0H+  +  H2  -» 

+ 

+ 

0 

CM 

X 

H 

1.5e-9 

[Ref.  4] 

H20+  +  e-  OH*  +  H 

6.5e-7  (300/Te)’/2 

[Ref.  4] 

0+  +  N2  N0+  +  N 

3.0e-10  (at  5eV) 

[Ref.  14] 

NO-*-  +  e-  -♦ 

N  +  0 

2.35e-8 

0-*-  +  NH3  -> 

0  +  NH3-*- 

1.2e-9 

[Ref.  15] 

NH3-*-  +  e-  ->  NH3 

1.0e-7 

(estimated) 

CIV  Reactions 

N2  N2-*- 

See  Table  3 

Photo-Ionization 

EUV  Wavelength 

Ionization  Cross 

Calculated 

[Angstrom] 

Section  [cm2] 

Reaction  Rate 

fcm3/sl 

H2  +  h-o  Ha-*-  +  e- 

800  -  630 

6e-18 

[Ref.  2] 

2.4e9 

1.8e-7 

630  -  460 

5e-18 

[Ref.  2] 

4.7e9 

1.5e-7 

460  -  370 

4e-18 

[Ref.  21 

0.63e9 

1 .2e-7 

N2  +  h-u  ->  N2-*-  +  e* 

800  -  630 

21e-18 

[Ref.  16] 

2.4e9 

6.3e-7 

630  -  460 

22e-18 

[Ref.  16] 

4.7e9 

6.6e-7 

460  -  370 

19e-18 

[Ref.  16] 

0.63e9 

5.7e-7 

370  -  270 

12e-18 

[Ref.  16] 

10.3e9 

3.6e-7 

270  -  205 

6e-18 

[Ref.  16] 

4.4e9 

1.8e-7 

NH3  +  hu  N 

H3-*-  +  e- 

1027-911 

18e-18 

[Ref.  17] 

11.61e9 

5.4e-7 

911  -800 

25e-18 

[Ref.  17] 

8.3e9 

7.5e-7 

800  -  630 

30e-18 

[Ref.  17] 

2.4e9 

9.0e-7 

630  -  550 

25e-18 

4.7e9 

7.5e-7 

Hgure  21  shows  the  computed  total  ion  distribution  along  the  solar  array.  The  array  is  located  5 
m  off  the  thnister  center  axis,  and  extends  about  10.S  m  downstream  of  the  thruster  exit  plane. 
For  this  simulation,  die  thnister  was  operating  into  the  wake  of  the  space  statitHi,  and  the 
simulatitxi  time,  which  corresponds  to  real  flow  dme  for  the  plume  exhaust,  was  3.5  ms.  The 
ambient  O'**  number  density  was  set  to  10^  cm*^  which  represents  a  high  value  for  the 
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ionosphere  in  LEO.  Total  ion  demines  along  the  array  can  be  seen  to  reach  values  as  high  as  3  x 
10^  cm'^.  which  is  a  13  times  increase  over  the  initial  O'*'  density.  Because  the  thruster  is  firing 
into  the  wake  at  3  km/s,  the  relative  velocity  of  the  plume  exhaust  with  respect  to  the  orbital  free 
stream  at  7.7  km/s  is  tmly  7.7  -  3  =  4.7  km/s;  this  is  to  slow  to  trigger  QV  effects.  The  reason  for 
the  local  plasma  build-up  is  that  the  relatively  slow  (~  3  km/s)  neutral  exhaust  plume  loses 
electrons  to  the  faster  (7.7  km/s)  O'*'  itms,  so  the  slow  moving  neutrals  become  a  slow  moving 
plasma.  Figure  22  shows  how  the  entering  from  the  left  has  been  depleted  downstream  of  the 
thruster  exit  plane  through  itm-molecule  reactitMis.  As  the  becomes  more  depleted 
(k)wnstream,  the  plasma  production  due  to  charge  exchange  decreases.  That  is  why  figure  21 
shows  peak  plasma  density  immediately  downstream  of  the  thruster  exit  plane. 

This  simulation  did  not  include  i^oto-ionizatitm,  even  though  the  numerical  code  includes  that 
option.  On  the  time  scales  we  are  looking  at-  (3.S  ms),  {dioto-itmization  is  too  slow  to  have  any 
effect  as  opposed  to  the  much  quicker  itm-molecule  reactions.  For  an  ambient  velocity  of  7.7 
km/s,  and  a  grid  cell  dimension  of  10  cm,  the  Courant  stability  ctmdition  requires  that  the  time 
step  be  less  than  10*^  s  for  the  numerical  scheme  to  be  stable.  For  this  simulation  we  used  a  time 
step  of  10*^  s.  If  the  photo-itmization  were  to  be  included,  the  velocity  in  the  Courant  stability 
condition  would  be  the  photon  velocity,  which  is  the  speed  of  light,  3x10^  km/s.  The  time  step 
must  then  be  on  the  orxier  10'^^  s,  and  the  number  of  ctxnputational  steps  to  simulate  3.5  ms 
increases  several  orders  of  magnitude.  To  keep  the  computational  effort  reasonable,  we  thus 
excluded  photo-ionization  from  all  the  3.5  ms  simulatitms. 

AtKNher  way  to  include  the  photo-ionization,  is  just  to  define  a  photo-ion  production  rate  for  each 
plume  spedes,  i.  e.,  a  fixed  number  of  neutral  molecules  {riioto-ionized  per  sectmd.  This  keeps 
the  photon  flux  out  of  the  calculation  all  togetter,  assuming  that  the  i^on  flux  is  the 
undepleted  solar  EUV  flux  everywhere  at  all  times.  Bernhardt  et  al.  [Ref.  3]  report  that  a  neutral 
barium  cloud  becomes  optically  thick  to  the  Sun's  ionizing  EUV  radiation  for  number  densities 
greater  than  Nmax  =  10^  cm'^-  This  effect  can  be  approximated  by  limiting  the  photo-ion 
production  rate  to  Pmax  =  (Ngiax)  (v)*  where  t)  is  the  ion  production  rate  per  sectmd.  This 
prevents  the  creation  of  unrealistically  large  plasma  densities  in  the  neutral  "pile-up"  region  along 
the  sdar  array  seen  in  figure  9.  A  neutral  effluent  flow  of  number  density  10^  m~3  becomes 
optically  tiiick  after  5  - 15  cm.  depending  on  the  neutral  spedes.  The  c^cal  depth  is  inversely 
proptmional  to  the  number  density.  Gearly,  the  optical  thickness  of  the  neutral  exhaust  must  be 
taken  into  account  to  compute  the  photo-ionization  prt^rly. 
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Figure  23  shows  the  total  i(»i  density  for  the  same  thruster  operating  into  the  ram.  The  ambiem 
O'*"  is  now  entering  from  the  right,  and  the  plume-ambient  relative  velocity  is  7.7  +  3  =  10.7 
km^.  This  velocity  is  greater  than  the  critical  itmization  velocity  for  N2,  so  we  can  expea  to  see 
some  CIV  effects  here.  The  maximum  itm  density  has  not  changed  noticeably  ffom  the  wake 
firing  case,  but  the  spatial  distribution  of  the  itxis  have.  The  maximum  charge  exchange  activity 
takes  place  10  m  downstream  of  the  thruster  where  the  O'*'  is  still  undepleted.  In  addition  there  is 
a  local  plasma  density  peak  10  m  downstream  of  the  thruster  and  about  3  m  off  the  thruster 
center  axis;  this  is  due  to  CIV  effects.  The  next  figure  24.  where  only  the  CIV  produced  N2'*'  ions 
are  plotted,  illuminates  this  even  better.  The  axial  velocity  of  the  neutrals  is  greatest  close  to  the 
center  axis,  (tecreasing  radially  towards  the  solar  array.  That  means  that  the  relative  velocity 
between  the  neutral  exhaust  and  the  ambient  O'**  is  greatest  close  to  the  center  axis,  which  again 
means  that  CIV  ionization  is  maximized  close  to  the  center  axis.  This  trend  is  opposed  by  the 
neutral  density  deperulence  of  the  CTV  process.  Because  OV  ionization  is  propotional  to  the 
neutral  density  (up  to  a  point),  the  density  dependent  part  of  the  CIV  rate  will  increase  with  the 
neutral  density  toward  the  solar  array  (see  figure  9  on  neutral  densities).  The  cmnbined  effea  of 
the  velocity  and  density  depetulent  parts  of  the  CIV  rate,  results  in  the  wedge  like  plasma  density 
peak  seen  in  figure  23. 

According  to  the  CTV  process  as  it  was  described  in  the  QV  section,  the  process  is  initiated  by 
fast  seed  ions  ("fast"  relative  to  the  ambient  plasma)  generated  by  some  imiization  process  such 
as  charge  exchange  between  die  fast  neutral  exhaust  and  the  slower  ambient  O'*'.  These  seed  ions 
go  on  to  heat  the  electrons  through  plasma  collective  effects,  and  the  high-energy  tail  of  the 
electrons  ionize  new  fast  neutrals.  There  is  presumably  an  amount  of  time  required  for  the 
necessary  number  of  seed  ions  to  be  produced,  before  the  CTV  process  really  takes  off.  In  our 
simulations,  we  have  assumed  that  the  production  of  fast  seed  ions  and  the  CIV  process  take 
place  simultaneously.  If  this  is  not  the  case  and  if  the  time  to  produce  sufficient  seed  ions 
exceeds  33  ms,  then  CTTV  ionization  will  not  happen  in  time  for  it  to  be  important  to  the  solar 
array  arcing  calculations. 

25  Ibf  Bipropellant  MMH-N2O4  Thruster 

We  repealed  the  ionization  calculations  fm-  the  bipropeUam  MMH-N204  thruster.  The 
bipropellant  thruster  plume  has  a  different  species  composition  than  the  moiK^rc^llant, 
consisting  of  H2.  N2,  H2O,  CO,  and  CO2-  The  biprop  plume  exhibits  the  same  set  of  ion- 
mtdecule  reactions  as  the  monoprop,  except  the  reactions  involving  NH3.  The  additional  ion- 
molecule  reacdons,  and  die  CTV  and  photo-ionization  reactions  experienced  by  the  bipn^  are 
listed  in  the  ftdlowing  taUe: 
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Table  5:  Ionization  Rates  Bipropeilant  Hydrazine 

Ion-Molecule  Reactions 

Reaction 

Reaction  Rate  fcm^/sl 

O-*-  +  H20  H 

H20'*'  +  0 

6.0e-10  (at  5  eV) 

[Ref.  14] 

H20  +  H20'‘'  HSO-*-  +  OH 

1.67e-9 

[Ref.  4] 

H30'''  +  e'  Neutrals 

6.3e-7  (300/Te)’/2 

[Ref.  4] 

O-*-  +  C02  H 

C02'‘'  +  0 

6.0e-10 

[Ref.  14] 

C02'*'  +  e-  ->  CO*  +  0* 

3.8e-7  (at  5  eV) 

[Ref.  3] 

O'*-  +  CO  -> 

CO'*'  +  0 

6.0e-1 1  (at  5  eV) 

[Ref.  14] 

02'''  +  e'  ^ 

0{1D)+0(1P) 

1.9e-7  (Te/300)’/2 

[Ref.  61 

CIV  Reactions 

C02  002-*- 

See  Table  3 

Photo-ionization 

EUV  Wavelength 

Ionization  Cross 

Calculated 

[Angstrom] 

Section  [cm^] 

Reaction  Rate 

fcm^/sl 

H20  +  h\)  h 

120-*-  +  e- 

1027-911 

14e-18  [Ref.  17] 

11.61e9 

4.2e-7 

91 1  -  800 

15e-18  [Ref.  171 

8.3e9 

4.5e-7 

CO  +  hu  CO'*'  +  e* 

885  -  800 

16e-18  [Ref.  16] 

6.4e9 

4.8e-7 

800  •  630 

20e-18  [Ref.  16] 

2.4e9 

6.0e-7 

630  -  460 

17e-18  [Ref.  16] 

4.7e9 

5.1e-7 

460  -  370 

14e-18  [Ref.  161 

0.63e9 

4.2e-7 

C02  +  hi)  C02'*'  +  e* 

885  -  800 

lisSLM 

6.4e9 

4.5e-7 

800  -  630 

2.4e9 

6.0e-7 

630  -  460 

4.7e9 

7.5e-7 

460  -  370 

0.63e9 

7.5e-7 

Figure  25  shows  the  total  ion  density  for  a  bipiopeilant  hydrazine  thruster  firing  into  the  wake. 
The  ambient  O'*'  number  density  is  still  10^  cm*^,  and  the  simulation  time  was  3.5  ms.  The  peak 
plasma  density  is  similar  to  the  monopropellant  run  (figure  21),  reaching  values  as  high  as  3  x 
10^  cm'^  along  the  solar  array  right  after  the  thruster  exit  plane.  Note  that  the  density  color 
scales  are  different  on  figures  21  and  25,  with  figure  25  displaying  better  resoluticHi.  The 
distinctive  plasma  density  peak  in  figure  25,  is  due  to  local  QV  effects.  The  neutral  axial  exhaust 
velocity  is  zero  in  this  region  (see  figure  13),  as  the  plume  is  deflected  by  the  solar  array  into  the 
backflow.  That  means  that  the  relative  velocity  between  the  neutrals  and  the  ambient  O'*'  is  just 
the  ambient  orbital  velocity  in  this  regitm,  or  7.7  km/s.  This  is  sufficient  to  CIV  iraiize  the  CO2 
present  in  the  biprop  plume,  which  has  a  critical  ionization  velocity  of  7.6  km/s.  The  CIV 
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icmization  of  the  CO2  ts  even  clearer  in  figure  26,  which  shows  the  number  density  of  the  C02'*' 
ions  (Mily. 

In  this  next  case,  we  simulate  the  25  Ibf  biprop  hydrazine  thruster  firing  into  the  ram  direction. 
Since  this  plume  contains  the  CIV  ionization  prone  CO2,  we  expect  to  see  substantial  CIV 
activity.  The  total  plasma  density  after  a  3.5  ms  simulation  is  plotted  in  figure  27,  with  the 
ambient  O'*'  entering  from  the  left.  The  plot  shows  an  enormous  plasma  density  increase,  with 
peak  densities  at  the  array,  10  m  downstream  of  the  thruster,  as  high  as  10^^  cm'^,  which  is  a 
four  order  of  magnitude  increase  over  the  ambient  O'*"  density  of  10^  cm'^!  For  this  case,  the 
relative  velocity  between  the  neutral  exhaust  and  the  free  stream  O'*'  exceeds  10  km/s  in  most  of 
the  simulated  region,  being  as  high  as  10.7  km/s  in  certain  areas.  Not  only  is  the  CO2,  with  a 
critical  velocity  of  7.6  km/s,  CIV  ionized  almost  everywhere,  but  N2,  with  a  critical  velocity  of 
10.3  km/s,  is  also  subject  to  CIV  ionizatitMi  in  some  places.  The  extent  of  the  CO2  ionization  is 
clearly  seen  in  figure  28,  where  the  CO2'*'  ion  density  is  plotted;  the  ion  population  is  almost 
exclusively  made  up  of  C02'*’. 

How  realistic  are  these  tremendous  CTV  iotuzatiem  gains?  In  the  following  some  of  the  numerical 
asstunptions  are  addressed,  and  comparisons  to  results  in  the  literature  are  made: 

•  As  we  have  already  mentioned,  we  do  not  include  any  initial  time  in  the  simulation  for  the 
production  of  fast  seed  ions  ("fast"  relative  to  the  ambient  plasma)  from  the  neutrals.  We 
assume  the  productitm  of  seed  ions  (from  charge  exchange)  and  the  CIV  process  to  start 
simultaneously.  If  this  seed  icm  production  time  is  substantial  compared  to  the  3.5  ms 
simulation  time,  the  QV  iortizatiem  will  not  kick  in  until  the  exhaust  effluents  have  blown 
past  the  10  m  solar  array,  thereby  lowering  the  plasma  density  at  the  array. 

•  For  the  runs  presented  here,  we  assumed  that  the  newly  QV  created  plasma  was  charge 
neutralized  immediately  by  electrons,  assuming  an  ample  supply  of  electrons  everywhere. 
This  in  it  itself  is  a  reasonable  assumptitxi,  however,  since  the  electron  density  is  part  of  the 
QV  ionization  rate  equation,  a  run-away  effea  is  created.  The  more  QV  ionized  plasma  that 
is  produced,  the  more  electrons  are  needed  to  charge  neutralize  the  plasma,  and  the  greater 
the  QV  ionization  rates  become.  As  the  equatiem  for  the  QV  ionization  rate  shows,  the  CIV 
plasma  increase  is  a  direct  function  of  the  electrrai  density  at  t=0,  and  exponentially 
proportional  with  time.  The  electron  density  at  t=0  was  tKX  held  ctmstant,  rather,  after  each 
time  step  of  the  numerical  simulation,  it  was  reset  to  the  current  electron  density.  After 
discussions  with  Biasca  [Ref.  13]  it  was  felt  thtu  this  was  a  reasonable  assumption. 
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•  The  CIV  ionizaticHi  did  have  an  upper  density  cutoff,  as  discussed  in  the  section  on  CIV,  and 
as  shown  in  the  numerical  rate  simulations  by  Biasca  et  al.  [Ref.  12], 

•  In  the  simularitm,  it  was  assumed  that  all  neutrals  that  exceeded  the  critical  ionization 
velocity,  were  subject  to  anomalous  CIV  ionization.  It  has  been  pointed  out  in  the  literature 
[Ref.  8]  that  the  critical  ionization  velocity  is  really  an  optimistic  minimum,  and  that  if  there 
is  no  external  power  source  to  provide  the  kit^tic  energy  to  the  neutrals,  then  the  CIV 
ionizaticm  rate  vrill  exhibit  a  less  than  unit  efficiency.  In  the  thruster  case  modeled  here, 
however,  the  neutral  exhaust  molecules  are  energized  by  the  thruster  activity,  and  we  feel 
that  the  CIV  ionization  efficiency  should  be  100%. 

•  Newell  [Ref.  8]  reports  that  in  the  Project  Porcupine  space  QV  experiment,  which  released 
barium  into  the  itmosphere  at  velocities  in  excess  of  its  critical  velocity  of  2.7  km/s,  the  rapid 
initial  ionizatitm  was  observed  to  be  30  ±  10%  of  the  total  neutral  population  that  had  a 
transverse  velocity  greater  than  2.7  km/s.  This  ionization  rate,  which  presumably  is  due  to 
av  effects,  is  truly  enormous!  How  does  that  compare  to  our  results?  For  the  biprop 
hydrate  thruster  firing  into  the  ram,  we  obtained  C02'*’  ion  densities  of  lO'®  cm*^.  The 
neutral  plume  exhaust  density  is  on  the  average  about  10^^  m*^,  and  if  the  CO2  molar 
fraction  is  0.06,  that  amounts  to  a  CO2  number  density  of  6  x  10^^  m*3  or  6  x  10^  i  cm*3.  If 
about  80%  of  the  CO2  molecules  have  velocities  in  excess  of  the  critical  velocity,  that  leaves 
us  with  about  5  x  10^  i  cm'^.  Let  us  furthermore  assume  that  the  QV  ionization  is  somewhat 
lower  than  that  observed  by  Porcupine,  or  20%.  The  CIV  itmization  then  becomes  20%  of  S 
X  lOi  i  cm'^,  or  10^  i  cm'^.  Amazingly  enough,  this  is  only  one  order  of  magnitude  greater 
than  the  itmizaticm  calculated  for  the  biptt^  thruster!  And  the  biprop  simulation  time  was 
(xily  3.S  ms,  whereas  the  Porcupine  release  lasted  several  sectmds.  The  good  correlation 
between  the  Porcupine  test  observations  aiKf  our  numerical  simulation  results  is  encouraging, 
and  indicates  that  CIV  effects  can  increase  the  local  plasma  density  at  space  system  surfaces 
several  orders  of  magnitude. 

•  Finally,  Biasca  et  al.  [Ref.  12]  have  calculated  plasma  densities  near  the  point  of  a  thruster 
firing  into  the  ram,  and  for  certain  CTV  reacticm  rates,  they  report  two  orders  of  magnitude 
local  plasma  density  increases.  Hence  the  simulated  QV  itmization  levels  presented  here,  fall 
between  the  Porcupine  observed  levels  arid  the  levels  calculated  by  Biasca  et  al.  [Ref.  12]. 
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Arcing  at  High-Voltage  Surfaces 

The  effluents  released  from  various  space  system  thmsters  generate  a  self-induced  neutral/plasma 
envircHiment  that  may  affect  spacecraft  operations.  This  study  focuses  on  the  interactions 
between  the  self-induced  envinmment  and  biased  high-voltage  surfaces  in  space.  As  an  example 
we  look  at  the  solar  arrays  that  are  baselined  for  the  Space  Station  Freedom.  The  75  kW  of  power 
needed  for  the  SSF  will  be  supplied  at  high  voltage  and  low  current  to  minimize  resistive  losses 
and  the  mass  of  cabling  and  harnesses.  The  arrays  will  be  operating  at  160  volt,  which  is 
significantly  higher  than  currently  used  arrays  with  voltage  drops  ranging  from  28  to  75  volt 
(Skylab). 

For  negatively  biased  solar  cell  interamnects,  it  is  observed  that  below  a  critical  voltage,  arc 
discharges  occur  on  the  solar  array.  Arcing  has  been  defined  as  a  sudden  current  pulse  much 
larger  than  the  ambient  current  coUection  and  typically  lasting  a  few  microseconds.  Negative 
voltage  thresholds  ranging  from  200  to  several  100  volt  have  been  observed  in  ground-based  tests 
[Ref.  18, 20].  The  flight  data  frcMn  the  plasma  interaction  experiments  FIX  I  and  FIX  II  show 
arcing  arcing  occurring  at  thresholds  as  low  as  -200  V  [Ref.  19,20]. 

In  a  microscopic  aicii^  mechanism  develqied  by  Hastings  et  al.  [Ref.  20]  and  Cho  and  Hastings 
[Ref.  21],  it  is  proposed  that  the  arcing  breakdown  occurs  by  itxiizatian  of  the  neutral  gas 
desorbed  frcnn  the  solar  array  dielectric  cover  glasses.  The  desorption  is  known  as  electron 
stimulated  desorption  (BSD)  with  the  electrons  being  emitted  from  the  solar  array 
interconiKCtors.  A  schematic  view  of  the  solar  array  process  leading  to  arcing  is  shown  in  figure 
29.  Note  that  only  the  dielectric  cover  glass  and  adhesive  are  shown;  the  solar  cell  is  left  out. 

This  is  because  the  solar  cell  is  a  semictmductor  with  a  voltage  drop  across  it  of  at  most  a  volt  or 
two,  wttich  is  negligible  in  the  charging  process.  The  proposed  arcing  mechanism  goes  as 
follows; 

•  Ambioit  positive  ions  are  attracted  to  the  -160  V  solar  arrays,  impinging  on  both  the 
interconnects  and  the  solar  cell  cover  glass.  The  ions  charge  the  dielectric  cover  glass 
positive  with  respect  to  the  negative  interccHinects,  leaving  the  cover  glass  sides  relatively 
uncharged.  A  strong  electric  field  is  thus  set  up  between  the  positive  cover  glass  surface  and 
the  interconnector,  the  field  is  maximized  at  the  triple  junction  point  where  the  conducting 
interconnector,  the  dielectric,  and  the  plasma  all  meet. 

•  The  electric  field  at  the  interconnector  surface  can  be  enhanced  by  whiskers  or  dielectric 
impurities;  the  enhancement  can  be  expressed  by  a  so-called  field  enhancement  factOT,  p. 


23 


Enhanced  field  emission  electrons  (EFEE)  from  the  interconncctor  impinge  on  the  cover 
glass  and  adhesive;  this  triggers  multiple  secondary  electron  emission  from  the  side  of  the 
dielectric  cover  glass  which  leaves  the  sides  positively  charged.  This  increases  the  electric 
field  across  the  triple  junction  point.  The  process  can  develop  very  rapidly  because  of  the 
strong  erqxHtential  dependence  of  the  EFEE  current  on  the  electric  field.  When  the  electric 
field  doubles,  the  emissitKi  currcnt  increases  by  ten  orders  of  magnitude. 

•  Neutral  gas  is  desorbed  into  the  gap  between  the  cells  due  to  heating  of  the  cover  glass  side 
by  the  electron  current.  The  EFEE  electrcm  current  flowing  through  the  neutral  cloud  can 
eventually  lead  to  a  discharge  like  a  surface  flash  over,  and  we  have  arcing. 

•  For  high  bias  voltage  (several  100  V  negative)  and  a  high  field  enhancement,  the  ion 
charging  time  governs  the  total  charging  time  and  hence  the  arc  rate.  The  ion  charging  time, 
and  therefore  the  arc  rate,  are  functions  of  the  local  ion  density.  For  low  voltage  and  low 
electric  field  enhancement,  the  total  charging  time  is  governed  by  the  EFEE  charging  time, 

«  The  arcing  voltage  threshold  depends  on  the  local  neutral  density,  but  is  independent  of  local 
itm  density. 


According  to  this  theory  [Ref.  20,21]  a  local  neutral  density  of  at  least  2.2  x  10^^  m*^  is  required 
to  trigger  the  Townsend  breakdown  and  arcing  at  the  array  shown  in  figure  29,  operating  at  -500 
V  with  respect  to  the  ambient.  Hence,  the  thruster  induced  neutral  densities  of  almost  10^^  m'^ 
that  we  calculated  along  the  array  in  figure  9  are  unlikely  to  trigger  arcing  at  the  SSF  -160  V 
arrays.  There  is  the  possibility  that  the  thruster  effluents  could  pile  up  in  the  gaps  between  the 
solar  cells,  and  that  this,  cmnbined  with  the  neutral  desorption  off  the  side  of  the  cover  glass 
could  lower  the  arcing  thresholds  to  voluges  less  negative  than  -SOO  V.  The  bottom  line, 
however,  is  that  SSF  solar  array  voltages  probably  are  too  low  to  induce  arcing. 

The  electric  power  requirements  for  future  space  systems  are  likely  to  increase,  and  this  will  most 
likely  be  accomplished  by  increasing  the  (grating  voltage  of  the  solar  arrays.  For  increasing 
voltage  bias,  it  becomes  important  to  predict  and  control  the  local  self-induced  neutral 
environments,  to  minimize  the  risk  of  arcing. 

If  array  voltages  and/or  neutral  densities  are  such  that  arcing  does  indeed  occur,  it  beccmies 
important  to  predict  the  arc  rates.  According  to  the  theory  by  Cho  and  Hastings  [Ref.  21  ],  the  arc 
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rate  is  exponentially  proportional  to  the  negative  voltage  fo*-  voltages  between  approximately  - 
SOO  and  -300  V.  Furthermore,  at  negative  voltages  below  approximately  -400  V,  the  arc  rate  is 
linearly  proportional  to  the  local  plasma  density.  Ihis  is  shown  in  figure  30  which  presents  the 
total  arc  rate  versus  the  bias  voltage  (both  tm  logarithmic  scales)  for  two  different  neutral 
dertsities.  An  array  containing  SOO  solar  cells  of  dimension  2  by  2  cm  (as  flown  on  FIX  II), 
operating  at  -SOO  V,  and  in  an  ambient  plasma  density  of  S  x  10^  m'^,  will,  according  to  the 
theory,  experience  an  overall  arc  rate  of  about  0.1  s*^.  If  the  local  plasma  density  is  increased 
two  orders  of  magnitude  to  S  x  10^  ^  m*^,  the  arc  rate  increases  correspondingly  to  10  s*^  Thus, 
the  four  order  of  magnitude  increase  in  local  plasma  density  computed  for  the  biprop  hydrazine 
thruster  firing  into  the  ram,  could  presumably  increase  the  arc  rates  overall  arc  rates  on  the  array 
to  an  astcHiishing  10^  s*^ !  Such  dire  possibilities  clearly  provide  an  incentive  to  analyze  the  self- 
induced  environment  thoroughly,  in  order  to  predict  show  stopping  arc  rates. 


CONCLUSION 

Various  thruster  activities  around  space  systems  in  LEO  generate  self-induced  neutral  and  plasma 
environments.  The  simulations  presented  in  this  report  indicate  that  for  certain  geometries  and 
thruster  conditions,  density  increases  along  the  solar  arrays  can  be  as  high  as  six  orders  of 
magnitude  for  the  neutrals,  and  four  orders  of  magnitude  for  the  local  plasma.  Hie  tremendous 
ionization  demonstrated  by  the  critical  ionization  velocity  mechanism  stressed  the  importance  of 
iiKluding  this  effect  when  total  ionization  of  the  neutral  effiuents  is  calculated.  On  the  time  rates 
of  interest  in  our  calculatitxis,  the  main  ionizatitm  mechanisms  were  itm-molecule  reactions 
(including  charge  exchange)  and  CTV;  photo-ionization  was  too  slow  to  be  important 

Despite  the  tremendous  self-induced  perturbaticxis  in  the  local  plasmaAieutral  distribution,  it  was 
concluded  that  the  density  increases  around  the  SSF  are  not  likely  to  trigger  arcing  at  the  160  V 
solar  arrays.  Arrays  currently  flown  on  other  space  systems  are  also  too  low  in  (grating  voltage 
to  exhibit  arcing.  Future  space  systems,  which  most  likely  will  be  (grating  arrays  at  increased 
voltage  biases,  will,  however,  have  to  worry  about  arcing.  Once  the  arcing  voltage  tiireshold  has 
been  reached,  it  becomes  important  to  predict  the  arcing  rates.  These  rates  are  a  strong  function 
of  the  local  plasma  density,  which  again  is  a  functicm  of  the  neutral  density.  In  order  to  predict 
arcing  thresholds  and  rates,  it  is  therefore  imperative  that  the  local  neutral  and  plasma 
environment  around  the  space  system  is  well  characterized. 
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APPENDIX  A 

The  Direct  Simulation  Monte  Carlo  Method 

The  direct  simulation  Monte  Carlo  (DSMQ  method  is  a  technique  for  computer  modeling  of  a 
real  gas  by  several  thousand  simulated  molecules.  It  is  not  feasible  to  track  each  individual 
molecule;  instead  the  method  tracks  simulated  molecules  where  each  simulated  molecule 
represents  as  many  as  10^^  real  molecules.  The  velocities,  temperatures,  and  positions  of  these 
molecules  are  stored  in  the  computer  and  are  modified  with  time  as  the  molecules  are  tracked 
through  representative  collisions  and  boundary  conditions  in  physical  space.  The  computational 
task  associated  with  the  direct  simulatirxi  becomes  feasible  when  the  gas  density  is  sufficiently 
low.  The  degree  of  rarefaction  of  a  gas  flow  has  traditic»ially  been  expressed  by  the  Knudsen 
number  which  is  the  ratio  of  the  mean  fiee  path  to  a  typical  dimension  of  the  flow  field.  A 
Knudsen  number  of  0. 1  has  been  quoted  as  the  boundary  between  ccmtinuum  and  transition  flow 
regimes.  An  improvement  in  determining  the  breakdown  of  continuum  flow  is  based  on  the  work 
by  G.  A.  Bird  [Ref.  22],  using  the  Bird  parameter.  P, 

P»U/pf  |dp/ds|, 

where  U  is  the  flow  velocity,  p  is  the  density,  f  is  the  collisitm  frequency,  and  s  is  the  streamline 
(fistance.  We  start  the  DSMC  calculatitnis  at  a  Bird  parameter  ctmtour  of  0.04.  Required  inputs  to 
the  simulation  are  flow  composition,  number  densities,  temperature,  and  vector  velocity.  The 
rarefied  flow  is  divided  into  grid  cells  dimensitxied  to  the  local  mean  free  path,  and  the  molecules 
are  marched  through  the  grid  over  time  steps  corresponding  to  the  inter-collisicm  time.  The  code 
is  made  computationally  efficient  by  uncoupling  fite  molecular  motion  and  collisions.  Within 
each  cell  all  molecules  are  considered  identical.  At  each  time  step  a  collisim  pair  is  picked  at 
random  and  the  collision  probability  is  calculated  based  on  cross  sectitm  and  relative  velocities. 
This  probability  is  compared  to  a  random  number  between  0  and  1  (hence  the  name  Montr 
Carlo),  and  if  the  probability  is  greater  than  the  number,  a  collision  occurs,  if  not,  a  new  collision 
pair  is  picked.  Enough  collisions  are  calculated  to  match  the  local  collision  frequency.  After  the 
collision(s),  the  mtriecules  are  advanced  according  to  their  relative  velocities,  and  new  position 
coordinates  are  computed. 

Ideally,  the  cell  dimensions  should  be  on  the  order  of  the  local  mean  free  path  or  less,  and  the 
time  step  should  be  on  the  order  of  the  time  between  collisions  or  less.  To  reduce  the  number  of 
cells  in  the  computational  grid,  however,  this  requirement  is  sometimes  relaxed  in  the  high 
density  flow  region  close  to  the  continuum-to-rarefied  transition  line.  Due  to  the  robustness  of 
the  DSMC  code,  this  assumption  has  no  noticeable  impact  m  the  flow  simulation  [Ref.  23]. 
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To  define  the  DSMC  computational  grid,  the  transition  line  at  P=0.04  is  divided  into  increments 
of  length  1  to  S  local  mean  free  paths.  The  increments  define  the  grid  cell  dimensions  at  the 
DSMC  stait  line.  As  the  grid  extoxls  into  the  rarefied  flow  region,  the  cell  dimensions  are  scaled 
roughly  with  the  local  mean  free  path.  Thus  for  tte  25  Ibf  monopropellant  hydrazine  thnister 
flowfield.  which  was  simulated  in  a  roughly  S  by  10  m  region,  we  defined  a  grid  of 
approximately  5,000  cells.  The  time  step  with  which  the  molecules  are  marched  is  roughly  equal 
to  the  inter-collision  time  at  the  start  line  (10*^  seconds).  The  time  step  is  kept  constant 
throughout  the  flowfield  so  several  steps  are  required  to  traverse  me  cell  in  the  rarefied  end  of 
the  flowfield. 

A  DSMC  simulation  does  not  converge  to  a  final  solution,  however,  the  simulation  becomes 
steady  and  statistically  more  correct  with  time.  Simulations  are  generally  run  until  each  grid  ceil 
has  been  sampled  about  10^  to  10^  times.  In  the  back  flow  of  a  thruster  plume  this  condition  has 
to  be  relaxed  at  least  an  order  of  magnitude.  For  the  25  Ibf  mono/biprop  hydrazine  thrusters 
modeled  here,  the  DSMC  simulations  tan  on  the  order  20  to  60  CPU  hours  m  our  Stardent  mini 
super-computer,  which  corresponds  to  several  100  milliseconds  of  real  flow  time.  The  DSMC 
provides  composition,  density,  vector  velocities,  temperatures  (translaiimal,  rotatimal, 
vibrational),  and  Mach  numbers  in  every  cell  of  the  flowfield.  By  manipulating  this  data, 
properties  such  as  pressure  and  mass  flux  can  be  calculated  anywhere  in  the  flowfield.  At 
surfaces  within  the  flow,  the  code  provides  pressure,  shear  stress,  and  incident  and  reflected 
energy  fluxes. 
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APPENDIX  B 

Flux  Solver  for  Computing  Local  Plasma  Density 

A  ccanputer  code  for  solving  a  set  of  ion  and  neutral  continuity  equations  was  developed  and  is 
enclosed.  The  code,  which  is  written  in  FORTRAN,  is  suitable  for  execution  on  a  PC.  An  initial 
neutral  distiibutitxi  data  file  (ARRAY.DAT)  containing  2-D  flowfield  grid  data  (number  density, 
vector  velocities)  is  required.  Currently,  the  total  neutral  density  and  the  average  velocities  are 
read.  The  fraction  of  various  species  in  the  i'^itial  neutral  distribution  are  set  in  the  code.  By 
sbghtly  changing  the  read  statement,  one  co^id  make  the  code  read  the  densities  and  velocities 
per  species,  since  this  information  is  available  frxxn  the  DSMC  simulation  anyway.  The  grid 
dimensions  are  currently  10  by  10  cm,  but  this  can  obviously  be  changed.  The  code  has 
provisitMis  for  both  monopropellant  and  bipropellant  hydrazine  thruster  plumes,  and  prompts  the 
user  for  ambient  plasma  density,  plume  direction  (firing  into  ram/wake),  photon  flux  direction, 
time  step,  and  simulation  time. 

The  code  includes  ion-molecule  ionization  (including  charge  exchange),  QV  ionization,  and 
photo-i(Hiization.  The  photo-ionizatiCHi  is  computed  based  cm  photon  flux  and  wavelength,  and 
the  ionization  rate  as  a  fimction  of  wavelength.  Due  to  the  friioton  speed,  the  Courant  stability 
condition  requires  that  the  numerical  time  step  be  less  than  10*'®  s.  At  such  a  short  time  step, 
extensive  run  times  are  required  to  simulate  flow  times  of  millisecond  duration.  In  afterthought, 
it  would  have  been  better  to  include  the  i^oto-ionization  as  an  option  tnily.  The  maximum  time 
step  without  photo-ionization  is  on  the  order  10*^  s.  Furthermore,  instead  of  computing  the 
photo-itxiization  yield  as  a  function  of  wavelength,  it  could  be  included  as  a  total  photo- 
ionization  production  rate,  e.  g.,  10'^  ions  produced  per  neutral  molecule  per  second.  This  is 
typically  done  in  other  codes. 

The  Crv  ionization  rates  are  derived  from  plots  provided  by  Biasca  [Ref.  13].  For  the  two 
species  where  CIV  ionization  rates  are  availidrle  (CO2  and  N2),  two  rates  which  are  valid  for  two 
different  tteutral  density  regimes  have  been  derived.  The  rates  are  of  course  valid  only  for 
plaana/neutral  relative  velocities  in  excess  of  the  critical  ionizaticxi  velocity.  The  maximum 
electron  current  available  to  charge  rKutralize  the  exponentially  produced  CIV  ions  is  limited  by 
the  availaMe  ambient  electron  flux. 

The  code  starts  with  the  initial  neutral  distribution,  and  then  turns  tm  the  various  ionization 
mechanisms.  After  each  time  step,  the  colunm  of  grid  cells  closest  to  the  thruster  get  their  neutral 
distribution  reset,  to  simulate  the  continuous  influx  of  neutrals  from  the  thruster.  Similarly,  the 
grid  ctriumn  on  the  side  of  the  computational  region  facing  the  ambient  O'**  flux  are  replenished 
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with  ambient  density  O'**  alter  each  time  step.  TIk  plasma  is  made  charge  neutral  alter  each  time 
step  by  setting  the  electron  density  equal  to  the  total  itm  density  (see  also  QV  limitations  in 
preceding  paragr^). 

The  final  i(m  densities  are  written  in  a  grid  format  to  an  ouq)ut  file,  lON.DAT. 
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Neutral  Particle  Beam  (NPB)  Platform  -  Artist's  Concept 
Note  Thruster  Plume  Exhaust  at  Rear  of  Vehicle 


FIGURE  1 
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Log  H2  Number  Density  [m-3] 
Chemical  Power  System  H2  Release 


FIGURE  3 


Log  H2  Static  Pressure  [torr] 
Chemical  Power  System  H2  Reiease 


FIGURE  4 


SDRC  I-DEAS  4.1:  System  Assembly 


SSF  Assembly  STAGE-2 


IDEAS**2  SSP  4.x  System  Assembly 
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FIGURE  8 


Logio  Constant  Number  Density  Contours  of  Neutral  Effluents 
1 00%  Thermal  Accommodation  at  Solar  Array 
25  Ibf  Monopropellant  Hydrazine  Thruster  HIGH  VOLTAGE 

Thrusting  past  solar  array  into  wake  SOLAR  ARRAY 
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FIGURE  9 


Constant  Absolute  Velocity  Contours  of  Neutral  Effluents  [m/sj 
25  Ibf  Monopropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  Into  wake 
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FIGURE  12 


Constant  Axial  Velocity  Contours  of  Neutral  Effluents  [m/s] 

Solar  Array  Temperature  =  333  K 

25  Ibf  Monopropellant  Hydrazine  Thruster  HIGH  VOLTAGE 

Thrusting  past  solar  array  Into  wake  SOLAR  ARRAY 


FIGURE  13 


nstant  Radial  Velocity  Contours  of  Neutral  Effluents  [m/s] 

Solar  Array  Temperature  =  333  K 

25  Ibf  Monopropellant  Hydrazine  Thruster  HIGH  VOLTAGE 

Thrusting  past  solar  array  Into  wake  SOLAR  ARRAY 


FIGURE  14 


Logio  Constant  Number  Density  Contours  of  Neutral  Effluents 

25  Ibf  Bipropellant  Hydrazine  Thruster  HIGH  VOLTAGE 

Thrusting  past  solar  array  into  wake  SOLAR  ARRAY 


FIGURE  15 


FIGURE  16 


Constant  Temperature  Contours  of  Neutral  Effluents  IK] 

2.5  Ibf  Monopropellant  Hydrazine  Thruster  H  l  GH  VOLT  AGE 

Thrusting  past  solar  array  Into  wake  SOLAR  ARRAY 
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RGURE 17 


FIGURE  18 


Constant  Temperature  Contours  of  Neutral  Effluents  [K] 

0%  Thermal  Accommodation  (Specular  Reflection) 

25  Ibf  Monopropellant  Hydrazine  Thruster  HI GH  VOLTAGE 

Thrusting  past  solar  array  Into  wake  SOLAR  ARRAY 


FIGURE  19 


Logio  Constant  Number  Density  Contours  of  Neutral  Effluents 
25  Ibf  Monopropellant  Hydrazine  Thruster 
0% Thermal  Accommodation  (Specular  Reflection)  HIGH  VOLTAGE 

Thrusting  past  solar  array  Into  wake  SOLAR  ARRAY 


FIGURE  20 


Logic  Constant  Number  Density  Contours  of  Electrons  [cm-Si 
(Electron  Density  =  Total  Ion  Density) 

25  Ibf  Monopropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  Into  wake 

III  3.5  millisecond  simulation  high  VOLTAGE 


JOm  FIGURE  21 


Log  10  Constant. Number  Density  Contours  of  0+  [cm-3] 

25  Ibf  Monopropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  Into  wake 

3.5  millisecond  simulation  HIGH  VOLTAGE 
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FIGURE  22 


Logio  Constant  Number  Density  Contours  of  Electrons  [cm-3] 

(Electron  Density  =  Total  Ion  Density) 

25  Ibf  Monopropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  into  ram 

3.5  millisecond  simulation  HIGH  VOLTAGE 
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FIGURE 23 


Logio  Constant  Number  Density  Contours  of  CO2+  [cm-3] 
25  Ibf  Bipropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  into  wake 
,,,  3.5  millisecond  simulation 


Logic  Constant  Number  Density  Contours  of  Electrons  [cm’3] 
Electron  Density  =  Total  Ion  Density 
25  Ibf  Bipropeltant  Hydrazine  Thruster 
Thrusting  past  solar  array  into  ram 
3.5  millisecond  simulation 
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FIGURE  27 


Logio  Constant  Number  Density  Contours  of  CO2*  [cm'3] 

25  Ibf  Bipropellant  Hydrazine  Thruster 
Thrusting  past  solar  array  into  ram 

3.5  millisecond  simulation  AGE 

II  SOLAR  ARRAY 
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FIGURE  28 


SCHEMATIC  VIEW  OF  SOLAR  ARRAY  CHARGING  PROCESS 

Enhanced  Field  Electron  Emission  (EFEE) 
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Log  Total  Arc  Rate  [a-l]  vs  Log  Negative  Bias  Voltage  [V] 
for  Two  Different  Densities  [Ref.  21] 


FIGURE  30 
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program  fluxes 

c - 

c  This  program  reads  a  2-D  array  of  neutral  densities  and  velocities 
c  (from  ARRAY.DAT),  and  computes  the  ionization  and  ion  densities 
c  over  time.  The  code  has  provisions  for  both  monoprop  hydrazine  and 
c  biprop  hydrazine.  The  user  is  prompted  for  ambient  0+  density  and 
c  direction,  direction  of  photon  flux,  time  step,  and  simulation  time, 

c  The  current  version  has  photo-ionization  included  only  for  the 

c  monoprop  plume;  on  short  time  scales  (milliseconds),  photo-ionization 
c  is  not  very  important  compared  to  the  other  ionization  mechanisms; 

c  ion-molecule  reactions  (including  charge  exchange)  and  Critical 

c  Ionization  Velocity  (CIV)  effects, 

c 

c  The  current  grid  has  cell  dimensions  10  by  10  cm,  which  for  orbital 

c  velocities  of  7-8  km/s  requires  time  steps  no  greater  than  le-6  s,  due  to 

c  the  Courant  stability  condition.  If  photo-ionization  is  included  (which 
c  it  is  for  the  monoprop  set-up,  the  time  step  can  be  no  greater 
c  than  le-10  s. 

c 

c  The  output  (ion  species  density  in  m-3)  is  written  to  the  file  ION. DAT. 
c 

c  fl(nx,ny)  :  x  fluxes 

c  gl(nx,ny)  :  y  fluxes 

c  nx,ny  :  grid  dimensions 

c  dt  :  stable  time-step 

c 

c - - - 

c  fphl(nx,ny)=flux  of  photons  (band  1)  in  x-direction  [ph/cm2-s) 

c  gphl(nx,ny)=flux  of  photons  (batid  1)  in  y-direction 

c  fph2(nx,ny)=flux  of  photons  (band  2)  in  x-direction 

c 
c 

c  fh2(nx,ny)  and  gh2(nx,ny)  =  h2-flux  in  x-  and  y-directions 

c  fn2(nx,ny)  and  gn2(nx,ny)  =  n2-flux 

c  fnh3(nx,ny)  and  gnh3(nx,ny)  =  nh3-flux 

c  fopl(nx,ny)  and  gopl(nx,ny)  =  0+-flux 

c  fohpl(nx,ny)  and  gohpl(nx,ny)  =  0H+-flux 

c  fh2opl(nx,ny)  and  gh2opl(nx,ny)  =  H20+-flux 

c  fnopl(nx,ny)  and  gnopl(nx,ny)  =  N0+-flux 

c  fnh3pl(nx,ny)  and  gnh3pl(nx,ny)  >  NH3-(— flux 

c  fn2pl(nx,ny)  and  gn2pl(nx,ny)  =  N2+-flux 

c  fh2pl(nx,ny)  and  gh2pl(nx,ny)  *  H2+-flux 

c  fo2pl(nx,ny)  and  go2pl(nx,ny)  =  02+-flux 

c  fcopl(nx,ny)  and  gcopl(nx,.ny)  =  C0+-flux 

c  fco2pl(nx,ny)  and  gco2pl(nx,ny)  =  C02+-flux 

c  fh3opl(nx,ny)  and  gh3opl(nx,ny)  =  H30+-flux 

c 

c  dh2(nx,ny)  s  H2  number  density  [#/cm3] 

c  dn2(nx,ny)  »  N2 

c  dnh3(nx,ny)  *  NH3 

c  dopl(nx,ny)  »  0+ 

c  dohpl(nx,ny)  »  0H+ 

c  dh2opl(nx,ny)  s  H20-*- 

c  dnopl(nx,ny)  =  N0+ 

c  dnh3pi(nx,ny)  =  NH3+ 

c  dn2pl(nx,ny)  =  N2+ 

c  dh2pl(nx,ny)  «  H2+ 

c  do2pl(nx,ny)  =  02+ 

c  dcopl(nx,ny)  *  C0+ 

c  dco2pl(nx,ny)  =«  C02+ 


dh3opl(nx,ny)  =  H30+ 

de(nx,ny)  =  Electrons 

dphl(nx,ny)  =  Photons  (Wave  band  1) 

dph2(nx,ny)  =  Photons  (Wave  band  2) 

dph7(nx,ny)  =  Photons  (Wave  band  7) 

parameter  (nx=90,ny=20)  ’MUST  BE  REPLACED  FOR  ALL  SUBROUTINES 

common/grid/dx,dy,dt  !IF  CHANGED!!! 

common/  f  luxphl  /  f  phi  (  nx ,  ny ) ,  gphl  ( nx ,  ny  ) 

common / f luxph2 / f  ph2 ( nx , ny ) , gph2 ( nx , ny ) 

common/ f luxph3/ f  ph3 ( nx , ny ) , gph3 ( nx , ny ) 

common / f luxph4 / f  ph4 ( nx , ny ) , gph4 ( nx , ny ) 

common/ f luxph5 / f  ph5 ( nx , ny ) , gphS ( nx , ny ) 

common/ fluxphb/ fph6 ( nx , ny ) , gph6 (nx , ny ) 

common/ fluxph7 / f ph7 ( nx , ny ) , gph7 (nx , ny ) 

common/ fluxH2/ fh2 ( nx , ny ) , gh2 ( nx , ny ) 

common/ fluxN2/fn2(nx,ny) ,gn2(nx,ny) 

common/ f luxNH3 / f nh3 ( nx , ny ) , gnh3 ( nx , ny ) 

common/ fluxH20/ fh2o ( nx , ny ) , gh2o(nx , ny ) 

common/ f luxCO /fco(nx,ny),gco(nx,ny) 

common/ f luxC02/ f co2 (nx, ny ) ,gco2(nx,ny) 

common/ f luxOpl/ f opl ( nx , ny ) , gopl (nx , ny ) 

common/ f luxOHpl/ f ohpl ( nx , ny ) , gohpl ( nx , ny ) 

common/ fluxH20pl/ fh2v.  '»1  (nx , ny ) , gh2opl(nx , ny ) 

common/fluxNOpl/fnopj  < rx,ny) ,gnopl(nx,ny) 

common/ f luxNH3pl / f nh3pi ( nx , ny ) , gnhipl ( nx , ny ) 

common/fluxN2pl/fn2pl(nx,ny) ,gn2pl(nx,ny) 

common/ fluxH2pl/fh2pl(nx,ny) ,gh2pl(nx,ny) 

common/flux02pl/fo2pl(nx,ny) ,go2pl(nx,ny) 

common/fluxCOpl/fcopl(nx,ny) ,gcopl(nx,ny) 

common/ fluxC02pl/fco2pl(nx,ny) ,gco2pl(nx,ny) 

common/ fluxH30pl/ fh3opl (nx , ny ) , gh3opl(nx , ny ) 

common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) , 

$  dopl ( nx , ny ) , dohpl ( nx , ny ) , dhiopl (nx , ny ) , dnopl ( nx , ny ) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opl(nx,ny) , 

$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloneu/ vxn( nx , ny ) , vyn ( nx , ny ) 
common / ve lo i on/ vxoxy ( nx , ny ) , vy oxy ( nx , ny ) 
common/ veloph/vxph( nx , ny ) , vyph ( nx , ny ) 
comraon/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/ rxratebi/cll,cl2,ci3,cl4,cl5,ci6,cl7,cl8,cl9 
common/phrateH2/cph3h2 , cph4h2 , cph5h2 
common/phrateN2/cph3n2 , cph4n2 , cph5n2 , cph6n2 , cph7n2 
common/ phrateNH3/ cphlnh3 ,  cph2nh3 ,  cph3nh3,  cph4nh3 
common/phrateH20/cphlh2o , cph2h2o 
common/phrateC0/cph2co, cph3co,cph4co,cph5co 
common/phrateC02/cph2co2 , cph3co2 , cph4co2 , cph5co2 
c ommon / s c r a t _a / s c_ax ( nx , ny ) , s c_ay ( nx , ny ) 
common/scrat_b/sc  bx(nx,ny) ,sc_by(nx,ny) 
common/scratn/sc_Exh2(nx,ny) ,sc_byh2(nx,ny) , 

$  sc_bxn2(nx,ny) ,sc_byn2(nx,ny) , 

$  sc_bxnh3(nx,ny) ,sc_bynh3(nx,ny) , 

$  sc_bxh2o(nx,ny) ,sc_byh2o(nx,ny), 

$  sc_bxco(nx,ny) ,sc_byco(nx,ny) , 

$  sc_bxco2 ( nx , ny ) , sc_by co2 ( nx , ny ) , 

$  sc  bxphl(nx,ny),sc_byphl(nx,ny),sc_bxph2(nx,ny), 


$  s  c_by ph2 ( nx , ny ) , s  c_bxph3 ( nx , ny ) , sc_byph3 ( nx , ny ) , 

$  s  c2bxph4 ( nx , ny ) , s  c”by ph4 ( nx , ny ) , s  c_bxph5 ( nx , ny ) , 

$  sc_by ph5  V  nx , ny ) , sc^bxphb ( nx , ny ) , sc_byph6 ( nx , ny ) , 

$  s  c”bxph7 ( nx , ny ) , s  c^by ph7 ( nx , ny ) 
common/civn2iiin/civn21 ,  civn22 ,  civn23 
common/civco21im/civco21 , civco22 , civco23 
common/civ/dciv(nx,ny) ,vcivn2(nx,ny) , vcivco2(nx,ny) 
common/ omega/ omegaN2 , omegaC02 

common/elrestl/elfluence, in2restrict,ico2restrict 
comraon/elrest2/dn2plcheck(nx,ny) ,dco2plcheck(nx,ny) 
common/dph/dphlamb , dph2amb , dph3amb , dph4amb , 

&  dph5amb,dph6amb,dph7amb 

dimension  xcoord(nx,ny) ,ycoord(nx,ny) ,dneutr(nx,ny) 
dimension  dph(nx,ny) 


open(unit=l,file='ion.dat' ) 
open(uni t=2 , f ile= ' array . dat ' ) 
do  10  i=l,nx 
do  10  j=l,ny 

read(2,*)xcoord(i, j ) ,ycoord(i, j) ,dneutr(i, j ) , 

&  vxn(i,j),vyn(i,j) 

dneutr(i, j)=dneutr(i, j)*l.e-6  iConvert  m-3  to  cm-3 
vxn(i, j)=vxn(i,j)*100.  IConvert  m/s  to  cm/s 

vyn(i, j)=vyn(i, j)*100. 

10  continue 

dx=10.  ![cm] 

dy=10.  ![cmj 
te=1500. 
veloelec=3.e7 
omegaN2=3429. 
omegaC02s2182. 
civn21=l.e6 
civn22=:1.5€7 
civn23=6.e7 
civco21=l.e6 
civco22=2.5e7 
civco23=7 .e7 

:  All  reaction  rates  in  cm3/s 

cl=1.7e-9 

c2=7.5e-8*(sqrt(300./te)) 
c3=1.5e-9 

c4=6.5e-7*(sqrt(300. /te)) 
c5=3.e-10 
c6=2.35e-8 
c7=1.2e-9 
c8=l.e-7 
c9*1090. 
cl0=545. 


'.electron  temp. 

.'Ambient  electron  velocity  in  LEO  [cm/s] 

!N2  gyro  frequency 
!C02  gyro  frequency 

! Density  limits  that  CIV  rates  are  valid  within 


!0+  +  H2  ->  OH-k  +  H 

!0H+  +  e-  ->  0*  +  H 

»0H+  +  H2  ->  H20+  +  H 

!H20+  +  e-  ->  OH*  +  H 

'.0+  +  N2  ->  N0+  +  N 

'.N0+  +  e-  ->  N  +  0(2.75  eV) 

»0+  +  NH3  ->  0  +  NH3+ 

.'NH3+  +  e-  ->  NH3  (just  estimated!.'!) 
!NH3/CIV 

!N2/CIV(estimated) 


The  remaining  rate  constants  are  for  biprop  mix  only; 


cll*6.e-10 
cl2»0. 

C13-2180. 

cl4=6.e-10 

cl5-3.8e-7 

cl6=1.9e-7*(sqrt(te/300. )) 

cl7=6.e-ll 

cl8=1.67e-9 

cl9=6. 3e-7*(sqrt(300. / te)) 


+  0  (at  5eV) 
CO  (at  5eV) 


!0+  +  H20  ->  H20+ 

!C02  +  0+  ->  02+  + 

!C02/CIV  (max) 

!0+  +  C02  ->  C02+  +  0  (at  5  eV) 
!C02+  +  e-  ->  CO*  +  0*  (at  5  eV) 
!02+  +  e-  ->  0(1D)+0(3P)  (Rees) 
!0+  +  CO  ->  C0+  +  0 
!H20  +  H20+  ->  H30+  +  OH 
!H30+  +  e-  ->  neutrals 


c  =  2.99793elO  !  Photon  speed  [cm/s] 

c  Photon  reaction  rates  are  (cross-section  (cm2]  *  c(cm/s]): 
c  H2  +  hv  ->  H2+  +  e- 

cph3h2=6.e-18  *  c  !800  -  630  Angstrom  (Bauer, Phys. Planet. Ion. ,p. 54) 

cph4h2=5.e-18  *  c  !630  -  460  A 

cph5h2=4.e-18  *  c  1460  -  370  A 

c  N2  +  hv  ->  N2+  +  e- 

cph3n2=21.e-18  *  c  !800  -  630  A  ("High  activity" ,Torr  et  al.,  1979) 

cph4n2=22.e-18  *  c  !630  -  460  A 

cph5n2=19.e-18  *  c  1460  -  370  A 

cph6n2=12.e-18  *  c  !370  -  270  A 

cph7n2=  6.e-18  *  c  !270  -  205  A 

c  NH3  +  hv  ->  NH3+  +  e- 

cphlnh3=18.e-18  *  c  !1027  -  911  A  (Really  TOTAL  absorption, Hudson, 71) 
cph2nh3=25.e-18  *  c  !911  -  800  A 

cph3nh3=30.e-18  *  c  !800  -  630  A 

cph4nh3=25.e-18  *  c  !630  -  550  (does  not  reach  460  A!) 

c  H20  +  hv  ->  H20+  +  e- 

cphlh2o=14.e-18  *  c  !1027  -  911  A  (Really  TOTAL  absorption, Hudson, 71) 
cph2h2o=15.e-18  *  c  !911  -  800  A 

c  CO  +  hv  ->  C0+  +  e- 

cph2co=16.e-18  *  c  !885  -  800  A  (does  not  reach  911  A) 

cph3co=20.e-18  *  c  !800  -  630  A  (Torr  et  al.,  1979) 

cph4co=17.e-18  *  c  !630  -  460  A 

cph5co=14.e-18  *  c  !460  -  370  A 

c  C02  +  hv  ->  C02+  +  e- 

cph2co2=15.e-18  *  c  !885  -  800  A  (does  not  reach  911  A) 

cph3co2=20.e-18  *  c  !800  -  630  A  (Torr  et  al.,  1979) 

cph4co2=25.e-18  *  c  !630  -  460  A 

cph5co2=25.e-18  *  c  1460  -  370  A 

c  Ambient  photon  density  for  the  different  wavelength  bands 
c  is  the  photon  flux  (#/cm2-s]  divided  by  c  [cm/s]: 
c  (Fluxes  are  from  S.J. Bauer,  "Physics  of  Planetary  Ionospheres" ,p. 47) 
dphlamb=11.61e9/c  ![#/cm3] 
dph2ambs8 . 3e9/c 
dph3amb=2.4e9/c 
dph4amb»4.7e9/c 
dph5amb*0.63e9/c 
dph6ambsl0. 3e9/c 
dph7ambs4 . 4e9/c 

write(*,*)'Monopropellant  (4831  H2,  29X  N2,  22X  NH3)  or' 
vrite(*,*)' biprop  fuel  (20X82,  31XN2,  32XH20,  lUCO,  6%C02):' 
write(*,*)' Enter  1  for  monoprop,  2  for  biprop:  ' 
read(*,*)  iprop 

write(*,*) 'Thruster  firing  into  ram  or  wake?' 
write(*,*)' Enter  1  for  wake,  -1  for  ram:  ' 
read(*,*)  ivox 

write(*,*) 'Enter  ambient  0+  density  (cm-3]:  ' 
read(*,*)  densoxyg 

write(*,*)'Enter  -1,  0,  or  1  for  negative,  no,  or  positive' 
write(*,*)' component  of  photon  velocity  X  component:  ' 
read(*,*)  ivxphoton 


vrite(*,*) 'Enter  -1,  0,  or  1  for  negative,  no,  or  positive' 
write(*,*)' component  of  photon  velocity  Y  component:  ' 
read(*,*)  ivyphoton 

rvox=f loat ( ivox) 
rvxphoton=float(ivxphoton) 
rvyphoton*float( ivyphoton) 

if  (iprop.eq.l)  then 
do  20  i=l,nx 
do  20  j=l,ny 

dh2(i, j )=0.48*dneutr(i,j) 
dn2(i, j )=0.29*dneutr(i, j ) 
dnh3( i , j )=0. 22*dneutr( i , j ) 
dh2o(i , j )=0. 
dco(i,j)=0. 
dco2(i,j)=0. 
dn2plcheck( i , j )=1 . 
dco2plcheck( i , j )=1 . 
dciv(i , j )=dneutr(i, j )/densoxyg 
vxoxy(i, j)=7.7e5*rvox 
vyoxy(i,j)=0. 
vxph(i,j)=c  *  rvxphoton 
vyph(i,j)=c  *  rvyphoton 
vci vn2 ( i , j ) =abs ( vxn ( i , j ) - vxoxy (i,j))/10.2e5 
20  continue 

else 

do  30  i«l,nx 
do  30  j=l,ny 

dh2(i,  j)=:0.20*dneutr(i,j) 
dn2(i, j )=0.31*dneutr(i, j ) 
dnh3(i,j)=0, 

dh2o( i , j )=0 . 32*dneutr(i , j ) 

dco(i , j )»0. ll*dneutr( i , j ) 

dco2(i,j)=0. 06*dneu  t  r ( i , j ) 

dn2plcheck( i , j )=1 . 

dco2plcheck(i, j )=1. 

dciv(i, j)*dneutr(i, j)/densoxyg 

rvqx= float(ivox) 

vxoxy(i, j )=7.7e5*rvox 

vyoxy(i,j)=0. 

vxph(i,j)*c  *  rvxphoton 

vyph(i,j)ac  *  rvyphoton 

vcivn2U.j)*abs(vxn(i, j)  -  vxoxy(i, j))/10.2e5 
vcivco2(i, j)sabs(vxn(i, j)  -  vxoxy (i,j))/7.7e5 
30  continue 
endif 

do  40  isl,nx  llnitial  conditions 
do  40  jal,ny 
dopl ( i , j )>densoxyg 
dohpl(i, j )*0. 
dh2opl(i,j)*0. 
dnopl(i,j)»0. 
dnh3pl(i,j).0. 
dn2pl(i,j).0. 
dh2pl(i,j).0. 
dco2pl(i, j)«0- 


do2pl(i, j)*0. 
dcopl(i, j)=0, 
dh3opl(i, j)=0. 
de(i, j)=densoxyg 
dphl ( i , j ) =dphlamb 
dph2 ( i , j ) =dph2amb 
d  ph3 ( i , j ) =d  ph3amb 
dph4(i,j)=d  phAamb 
dph5( i , j )=dph5amb 
dph6 ( i , j ) =dph6amb 
d  ph7 ( i , j ) =dph7 amb 
AO  continue 
step=0. 
tiine=0. 

write(*,*) 'Enter  time  step  [seconds]:  ' 
read(*,*)  dt 

write(*,*) 'Enter  simulation  time  [#  of  time  steps];  ' 
read(*i*)isimtime 

The  max  ambient  electron  fluence  available  per  time  step: 
(electron  density)*(electron  velocity)*( time) 
elfluence  *  densoxyg*veloelec*dt 

if  (iprop.eq.2)  go  to  100 
do  50  n=l,isimtime 
call  neutral 
call  oplus 
call  ohplus 
call  hloplus 
call  noplus 
call  nhiplus 
call  n2pius 
call  h2plus 
call  photon 

do  170  i=2,nx-l 
do  170  j=2,ny 

dneutr(i, j)*dh2(i,3)  +  dn2(i,j)  +  dnh3(i,j) 

&  +dh2o(i,j)  +  dco(i,j)  +  dco2(i,j) 

dph(i, j)=dphl(i, j)  +  dph2(i,j)  +  dph3(i,j)  +  dphA(i,j) 

&  +  dph5(i,j)  +  dph6(i,j)  +  dph7(i,j) 

if  (iprop.eq.l)  then 

write(3,*)xcoord(i, j),ycoord(i, j)»dneutr(i, j),de(i, j), 

&  dopK i , j ) , dohpl( i , j ) , dh2opl( i , j ) ,dnopl( i , j ) , dnh3pl( i , j ) , 
&  dn2pl(i,j),dh2pl(i,j),dph(i,j) 
else 

write(3,*)xcoord(i, j),ycoord(i, j),dneutr{i,j),de(i, j), 

&  dopl(i,j),dohpl(i,j),dh2opl(i,j),dnopl(i,j)»dco2pl(i,j), 
&  dn2pl(i,j),do2pl(i,j),dcopl(i,j),dh3opl(i,j),dph(i,j) 
end  if 

170  continue 

do  60  j»l,ny 
if  (ivox.eq.l)  then 
dopl(l, j )*densoxyg 
else 

dopl ( nx , j ) >densoxyg 
end  if 

if  (ivxphoton.eq.l)  then 
dphl ( 1 , j ) «dphlamb 
dph2(l, j )»dph2amb 


dph3 ( 1 , j ) =dph3amb 
dph4( 1 , j )=dph4amb 
dph5 ( 1 , j ) =dph5amb 
dph6(l, j )=dph6amb 
dph7(l,  j  )=:dph7amb 
end  if 

if  (ivxphoton.eq.-l)  then 
dphl ( nx , j ) =dphlamb 
dph2(nx, j )=dph2amb 
dph3(nx, j )=dph3amb 
dph4(nx, j )=dph4amb 
dph5.(nx,  j  )=dph5amb 
d  ph6 ( nx , j ) =d  phbamb 
dph7 ( nx , j ) =dph7amb 
end  if 

dh2(l, j)=0.48*dneutr(l, j ) 
dn2(l, j )=0. 29*dneutr(l, j ) 
dnh3(i, j )=0.22*dneutr(l, j ) 
dh2o(l,j)=0. 
dco(l, j)*0. 
dco2(l,j)=0. 

dciv(l, j)=dneutr(l, j )/densoxyg 
60  continue 

do  70  i  =  l,nx 
if  (ivyphoton.eq.l)  then 
dphl ( i , 1 )=dphlamb 
dph2 ( i , 1 ) »dph2amb 
dph3(i , l)*dph3amb 
dph4 ( i , 1 ) adph4amb 
dphS  (  i ,  1 )  «:dph5aiBb 
dph6( i , 1 )sdph6amb 
dph7(i ,l)=dph7amb 
end  if 

if  (ivyphoton.eq.-l)  then 
dphl(i ,ny}=dphlamb 
dph2(i,ny)-dph2amb 
dph3  (  i ,  ny  )  >dph3ainb 
dph4( i , ny)=dph4amb 
dph5( i , ny )=dph5amb 
dph6 ( i , ny ) adphbamb 
dph7 ( i , ny )=dph7amb 
end  if 

70  continue 

do  80  is2,nx 

dh2(i , l)»0.48*dneutr(i , 1) 

dn2(i,l)*0.29*dneutr(i,l) 

dnh3 ( i , 1 ) «0 . 22*dneu  t  r ( i , 1 ) 

dh2o(l,j)*0. 

dco(l,j)=0. 

dco2(l, 3)=0. 

dciv(i, l)adneutr(i, l)/densoxyg 
80  continue 

do  90  isl,nx 
do  90  jal,ny 

if  (dopl(i,j).lt.0)  dopl(i,j)*0 

i » j  )=dopl( i ,  j  )+dohpl(i ,  j  )+dh2opl(i ,  j  )+dnopl(i ,  j ) 


&  +dnh3pl(i,j)+dn2pl(i,j)+dh2pl(i,i) 
90  continue 

time=time+dt 
step=step+l . 

50  continue 
go  to  155 
100  continue 

do  110  n=l,isimtime 
call  neutralbi 
call  oplusbi 
call  ohplus 
call  h2oplusbi 
call  noplus 
call  co2plusbi 
call  n2plus 
call  o2plusbi 
call  coplusbi 
call  hSoplusbi 
call  photon 

do  120  j=l,ny 
if  (ivox.eq.l)  then 
dopl ( 1 , j ) =densoxyg 
else 

dopl(nx, j )=densoxyg 
endif 

if  (ivxphoton.eq.l)  then 
dphl(l, j )»dphlamb 
dph2  ( 1 ,  j  )  »dph2ainb 
dph3 ( 1 , j ) »dph3amb 
dph4( 1 , j )*dph4amb 
dph5(l, j)=dph5arab 
dph6(l, j )=dph6amb 
dph? ( 1 , j ) =dph7amb 
endif 

if  (ivxphoton.eq.-l)  then 
dphl(nx ,  j  )x>dphlamb 
dph2(nx,j)=dph2amb 
dph3(nx, j )»dph3amb 
dph4(nx,  j  )>:dph4ainb 
dph5(nx,j).dph5amb 
dph6 ( nx , j ) «d  ph6amb 
dph7  ( nx ,  j  )  =:dph7amb 
endif 


dh2(i , l)«0.20*dneutr(i , 1) 
dn2(i , l)a0. 31*dneutr(i , 1) 
dnh3(i,l)-0. 

dh2o(i,l)»0.32*dneutr(i,l) 

dco(i,l)-0.11*dneutr(i,l) 

dco2(i , l)*0.06*dneutr(i , 1) 
dciv(l , j )*dneutr(l , j )/densoxyg 
120  continue 

do  130  i  .  l,nx 
if  (ivyphoton.eq.l)  then 
dphl ( i , 1 )«dphlamb 
dph2 ( i , 1 )«dph2afflb 
dph3(i,l)sdph3amb 


dphA  (  i ,  1 )  =dph4ainb 
dph5 ( i , 1 ) sdphSamb 
dph6(i , l)=dph6amb 
dph7 ( i , 1 ) =dph7  amb 
end  if 

if  (ivyphoton.eq. -1)  then 
dphl ( i , ny )=dphlamb 
dph2(i ,ny)=dph2amb 
dph3( i , ny)=dph3amb 
dphA  ( i ,  ny  )  ^^dphAamb 
dph5(i ,ny)=dph5amb 
dph6(i,ny)=dph6amb 
dph7 ( i , ny )=dph7amb 
end  if 

130  continue 

do  lAO  i=2,nx 
dh2(i , l)=0.20*dneutr(i, 1) 
dn2(i , l)=0.31*dneutr(i, 1) 
dnh3(i,l)=0. 

dh2o(if l)=0.32*dneutr(i,l) 
dco(i,i)=0.11*dneutr(i,i) 
dco2(i, i)=0.06*dneutr(i , 1) 
dci v( i , 1 )=dneutr ( i , 1 ) /densoxyg 

lAO  continue 

do  150  i=l,nx 
do  150  j=l,ny 

if  (dopl(i,j).lt.0)  dopl(i,j)=0 

de(i, j )adopl(i, j)+dohpl(i, j)+dh2opl(i» j)+dnopl(i, j) 

&  +dn2pl(i , j )+do2pi(i , j )+dcopi(i , j )+dco2pi(i , j ) 

&  +dh3opl(i,j) 

150  continue 

time=time+dt 

step=step+l. 

110  continue 

155  continue 

do  160  i=2,nx-l 
do  160  j=2,ny 

dneutr(i, j)*dh2(i,j)  +  dn2(i,j)  +  dnh3(i,j) 

&  +dh2o(i,j)  +  dco(i,j)  +  dco2(i,j) 
dph(i,j)=dphl(i,j)  +  dph2(i,j)  +  dph3(i,j)  +  dphA(i,j) 

&  +  dph5(i,j)  +  dph6(i,j)  +  dph7(i,j) 

if  (ipcop.eq.l)  then 

write(l,*)xcoord(i,j),ycoord(i, j),dneutr(i, j),de(i,j), 

&  dopl ( i , j ) , dohpl ( i , j ) , dh2opl ( i , j ) , dnopl ( i , j ) , dnh3pl ( i , j ) , 
&  dn2pl(i,j),dh2pl(i,j),dph(i,j) 

else 

wr i te( 1 , * )xcoord ( i , j ) , y coord ( i , j ) , dneu t  r (i,j),de(i,j), 

&  dopl(i,j),dohpl(i, j),dh2opl(i,j),dnopl(i, j),dco2pl(i,j)» 
&  dn2pl(i,j),do2pl(i, j).dcopl(i,j),dh3opl(i,j)»dph(i,j) 
end  if 

160  continue 
end 


subroutine  photon 


parameter  (nx=90,ny=20) 
common/gr id/dx , dy , d t 

common/ fluxphl / f phi ( nx , ny ) , gphl ( nx , ny ) 
common / f luxph2 / f  ph2 ( nx , ny ) , gph2 ( nx , ny ) 
common/ fluxphS/ fph3 ( nx , ny ) , gph3 ( nx , ny ) 
common/ fluxphA/ f ph4 ( nx , ny ) , gph4 ( nx , ny ) 
common/ fluxphS/ fph5 ( nx , ny ) , gph5 ( nx , ny ) 
common/  f  luxph6  /  f  ph6  (  nx ,  ny ) ,  gph6  ( nx ,  ny  ) 
common / f luxph? / £  ph7 ( nx , ny ) , gph? ( nx , ny ) 
common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) , 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opi(nx,ny),dnopl(nx,ny) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny),dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny),dh3opi(nx,ny) , 
$  dphl ( nx , ny ) , dph2 ( nx , ny ) , dph3 ( nx , ny ) , dph4 ( nx , ny ) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloph/vxph(nx,ny) ,vyph(nx,ny) 
common/ phrateH2/cph3h2 , cph4h2 , cph5h2 
common/ phra  teN2 / cph3n2 , cph4n2 , cph5n2 , cph6n2 , cph7n2 
common/phrateNH3/cphlnh3 , cph2nh3 , cph3nh3 , cph4nh3 
common/phrateH20/cphlh2o,cph2h2o 
common/phrateC0/cph2co , cph3co , cph4co, cph5co 
common/ phrateC02/cph2co2 , cph3co2 , cph4co2 , cph5co2 
common/ s  c  ra  t_a/ s  c_ax ( nx , ny ) , sc_ay ( nx , ny ) 
common/scrat_b/sc  bx(nx,ny) ,sc_by(nx,ny) 
common/scratn/sc_Exh2(nx,ny) ,sc_byh2(nx,ny) , 

$  sc_bxn2(nx,ny) ,sc_byn2(nx,ny) , 

$  sc”bxnh3(nx,ny) ,sc_bynh3(nx,ny) , 

$  s  c~bxh2o ( nx , ny ) , s  c~by h2o ( nx , ny ) , 

$  sc”bxco(nx,ny),sc_Eyco(nx,ny), 

$  sc~bxco2 ( nx , ny ) , sc_byco2 ( nx , ny ) , 

$  sc_bxphl(nx,ny) ,sc_byphl(nx,ny) ,sc_bxph2(nx,ny) , 

$  sc”byph2(nx,ny) ,sc_bxph3(nx,ny),sc_byph3(nx,ny) , 

$  s  c_bxph4 ( nx , ny ) , s  c_by ph4 ( nx , ny ) , sc  bxphS ( nx , ny ) , 

$  s  c_byph5 ( nx , ny ) , s  c~bxph6 ( nx , ny ) , sc_byph6 ( nx , ny ) , 

$  sc_bxph7(nx,ny) ,sc  byph7(nx,ny) 
common /dph/dph lamb , 3ph2amb , dph3amb , dph4amb , 

&  dph5amb,dph6amb,dph7amb 

initialize  scratch  arrays 

do  100  i=:l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  =  0,0 
sc_bxphl(i, j)  =  0.0 
sc_byphl(i , j )  =  0.0 
sc_bxph2(i,j)  =  0.0 
sc_byph2(i, j)  «  0.0 
sc_bxph3(i , j )  »  0.0 
sc_byph3(i , j )  =  0.0 
sc_bxph4(i,j)  *  0.0 
sc_byph4(i, j )  =  0.0 
sc_bxph5(i, j)  =0.0 
sc_byph5(i,j)  =  0.0 
sc_bxph6(i, j)  =0,0 
sc_byph6(i , j )  =  0.0 
sc_bxph7(i , j )  =  0.0 
sc  byph7(i,j)  =  0.0 


continue 

nxl=nx-l 

nyl*ny-l 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxph(i+l,j)  +  vxph(i,j)) 
continue 
do  lAO  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyph(i,j+l)  +  vyph(i,j)) 
continue 

do  150  i=l,nxl 
do  150  j=l,ny 

if(sc_ax(i,j).ge.0.)then 
sc_bxphl(i,j)=dphl(i,j) 
sc_bxph2(i,j)=dph2(i,j) 
sc_bxph3( i , j )-dph3 ( i , j ) 
scbxphA  < i , j ) =dph4( i , j ) 
sc_bxph5(i,j)=dph5(i,j) 
sc_bxph6(i,j)=dph6(i,j) 
sc_bxph7(i,j)-dph7(i,j) 
else 

sc_bxphl(i,j)*dphl(i+l,j) 
sc_bxph2 ( i , j )=dph2 ( i+1 , j ) 
sc_bxph3(i,j).dph3(i+l,j) 
sc_bxph4(i,j).dph4(i+l,j) 
sc_bxph5( i , j )«dph5( i+1 , j ) 
sc_bxph6(i,j).dph6(i+l,j) 
sc  bxph7(i,j).dph7(i+l,j) 
end  if 
continue 


do  160  i=l,nx 
do  160  j»l,nyl 
if(sc_ay(i,j).ge 
sc_byphl(i,j) 
sc_byph2(i,j) 
sc_byph3(i,j) 
sc_byph4(i,j) 

sc_byph5(i,j) 

sc_byph6(i,j) 

sc_byph7(i,j) 

else 

sc_byphl(i,j) 
sc_byph2(i,j) 
sc_byph3(i,j) 
sc_byph4(i,j) 
sc_byph5(i,j) 
sc_byph6(i,j) 
sc_byph7(i,j) 
end  if 
continue 


.0. )then 
=  dphl(i,j) 

=  dph2(i,j) 

*  dph3(i,j) 

=  dph4(i,j) 

=  dph5(i,j) 

=  dph6(i,j) 

=  dph7(i,j) 

-  dphl(i,j+l) 
=  dph2(i,j+l) 

*  dph3(i,j+l) 
=  dph4(i,j+l) 
«  dph5(i,j+l) 
»  dph6(i,j+l) 

*  dph7(i,j+l) 


do  170  isl,nxl 
do  170  j*i,ny 

fphl(i,j)  .  sc_ax(i,j)  *  sc_bxphl(i, j) 
fph2(i,j)  .  sc_ax(i,j)  *  sc_bxph2(i,j) 


fph3(i,j)  =  sc_ax(i,j)  *  sc_bxph3(i, j) 

fph4U»j)  =  sc_axU*j)  *  sc_bxph4(i,j) 

fph5(i,j)  =  sc_ax(i,j)  *  sc_bxph5(i , j ) 

fph6(i,j)  =  sc_ax(i,j)  *  sc_bxph6(i, j ) 

fph7(i,j)  =  sc_axU,j)  *  sc_bxph7(i , j ) 

170  continue 
c 

do  180  i*l,nx 
c 

do  180  j*l,nyl 

gphl(i,j)  =  sc_ay(i,j)  *  sc_byphl(i, j) 

gph2(i,j)  =  sc_ay(i,j)  *  sc_byph2(i, j ) 

gph3(i,j)  =  sc_ay(i,j)  *  sc_byph3(i,j) 

gph4(i,j)  =  sc_ay(i,j)  *  sc_byph4(i, j ) 

gph5(i,:j)  =  sc_ay(i,j)  *  sc_byph5(i,j) 

gph6(i,j)  =  sc_ay(i,j)  *  sc_byph6(i,j) 

gph7(i,j)  =  sc_ay(i,j)  *  sc_byph7(i, j) 

180  continue 
c 

do  200  i=2,nx 

do  200  j*2,ny 

dphl(i,j)  =  dphl(i,j)  -dt*(fphl(i,j)-£phl(i-l, j))>'dx 
+  -dt*(gphl(l,j)-gphl(i,j-l))/dy 

+  -dt*dphl(i,j)*(cphlnh3*dnh3(i,j)  +  cphlh2o*dh2o(i, j )) 
c  To  prevent  numerical  overflow: 

if  (dphl(i,j).lt,l.e-12)  dphl(i,j)=0. 

dph2(i,j)  =  dph2(i,j)  -dt*(fph2(i,j)-fph2(i-l, j))/dx 
^  -dt*(gph2(l,j)-gph2(i,j-l))/dy 

+  -dt*dph2(i, j)*(cph2ah3*dnh3(i, j)  +  cph2h2o*dh2o(i, j ) 

+  +  (cph2co*dco(if j )  +  cph2co2*dco2(i,3))*(85./110. )) 

if  (dph2(i,j).lt.l.e-12)  dph2(i,j).0. 

c  Photo-ionization  yield  of  CO  and  C02  within  wavelength  band  2 
c  (911  -  800  Angstrom)  is  reduced  since  ionization  cut-off 
c  occurs  at  885  A.  (911-800)/(885-800)-(85./lll. ) 

dph3(i,j)  =  dph3(i,j)  -dt*(fph3(i,j)-fph3(i-l,j))/dx 
+  -dt*(gph3(i,j)-gph3(i,j-l))/dy 

+  -dt*dph3(i,j)*(cph3h2*dh2(i,j)  +cph3n2*dn2(i, j ) 

+  +cph3nh3*dnh3(i, j ) 

+  +  cph3co*dco(i,j)  +  cph3co2*dco2(i , j ) ) 

if  (dph3(i,j).lt.l.e-12)  dph3(i,j)-0. 

dph4(i,j)  .  dph4(i,j)  -dt*(fph4(i,j)-fph4(i-l,j))/dx 
+  -dt*(gph4(i,j)-gph4(i,j-l))/dy 

+  -dt*dph4(i,j)*(cph4h2*dh2(i,j)  +cph4n2*dn2(i, j ) 

+  +cph4nh3*dnh3(i> j )*(80./170. ) 

-t-  +  cph4co*dco(i,  j)  +  cph4co2*dco2(i,  j)) 

if  (dph4(i,j).lt.l.e-12)  dph4(i,j)-0. 
c  Note  that  photo-ionization  yield  of  NH3  within  wavelength  band  4 
c  is  reduced  by  a  factor  of  80/170.  This  is  because  550  Angstrom  is 

c  the  cut-off  wavelength  for  ionization,  whereas  the  band  width 

c  is  630  -  460  A.  (630-550)/(630-460)  -  80/170. 

dph5(i,j)  .  dph5(i,j)  -dt*(fph5(i,j)-fph5(i-l, j))/dx 
+  -dt*(gph5(i,j)-gph5(i, j-l))/dy 

+  -dt*dph5(i,j)*(cph5h2*dh2(i,j)  +cph5n2*dn2(i, j ) 

+  +  cph5co*dco(i, j)  +  cph5co2*dco2(i, j )) 


if  (dph5(i,j).lt.l.e-12)  dph5(i,j)«0. 


dph6(i,j)  =  dph6(i,j)  -dt*(fph6(i,j)-fph6(i-l,j))/dx 
+  -dt*(gph6(i,j)-gph6(i, j-l))/dy 

+  -dt*dph6(i, j )*(cph6n2*dn2(i, j ) ) 

if  (dph6(i,j).lt.l.e-12)  dph6(i,j)=0. 

dph7(i,j)  =  dph7(i,j)  -dt*(fph7(i,j)-fph7(i-l,j))/dx 
+  -dt*(gph7(i,j)-gph7(i,j-l))/dy 

+  -d t*dph7 ( i , j )* ( cph7n2*dn2( i , j ) ) 

if  (dph7(i,j). It. 1.6-12)  dph7(i,j)=0. 

200  continue 
return 
end 
c 
c 

subroutine  neutral 
c 

parameter  (nx=90,ny=20) 

comfflon/grid/dx,dy,dt 

common/  f  luxH2  /  f  h2  (nx ,  ny  ) ,  gh2  (nx ,  ny  ) 

common/ fluxN2/ fn2 ( nx , ny ) , gn2 ( nx , ny ) 

common/ f  luxNH3  /  fnh3  (  nx ,  ny  ),  gnh3  (  nx ,  ny) 

common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) f 
$  dopl ( nx , ny ) , dohpl ( nx , ny ) , dhiopl ( nx , ny ) , dnopl ( nx , ny ) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opl(nx,ny), 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dphA(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloneu/ vxn ( nx, ny ) ,vyn(nx,ny) 
common/ veloion/vxoxy(nx , ny ) , vyoxy (nx , ny ) 
common/rxra  te/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/phrateH2/cph3h2 , cphAhZ , cph5h2 
common/phrateN2/cph3n2 , cph4n2 , cph5n2 , cph6n2 , cph7n2 
common/phrateNH3/cphlnh3, cph2nh3,cph3nh3, cph4nh3 
common/ s era t_a/ sc_ax ( nx , ny ), s c_ay( nx , ny ) 
common/scrat_b/sc  bx(nx,ny) ,sc_by(nx,ny) 
common/scratn/sc_Exh2(nx,ny) ,sc_byh2(nx,ny) , 

$  sc_bxn2(nx,ny) ,sc_byn2(nx,ny) , 

$  sc_bxnh3(nx,ny) ,sc_bynh3(nx,ny) , 

$  s  c_bxh2  0 ( nx , ny ) , s c~byh2o ( nx , ny ) , 

$  sc_bxco(nx,ny),sc_Byco(nx,ny), 

$  sc_bxco2(nx,ny) ,sc_byco2(nx,ny), 

$  s c_bxphl ( nx , ny ) , s c~by phi ( nx , ny ) , sc_bxph2 ( nx , ny ) , 

$  sc_byph2 ( nx , ny ) , sc_bxph3 ( nx , ny ) , sc”byph3(nx , ny ) , 

$  sc_bxph4(nx,ny) ,sc_byph4(nx,ny),sc”bxph5(nx,ny) , 

$  sc_by ph5 ( nx , ny ) , s  c_bxph6 ( nx , ny ) , sc^byphb ( nx , ny ) , 

$  sc_bxph7(nx,ny) ,sc_byph7(nx,ny)  ” 
common/civn2iim/civn21 , civn22 , clvn23 
common/civ/dciv(nx,ny) , vcivn2(nx,ny) ,vcivco2(nx,ny) 
comfflon/omega/omegaN2 , omegaC02 

initialize  scratch  arrays 

do  100  i>l,nx 
do  100  j»l,ny 
sc_ax(i,j)  =  0.0 
sc  ay(i,j)  «  0.0 


c 

c 


sc_bxh2(i,j)  =  0.0 
sc_byh2(i,j)  =0.0 
sc_bxn2(i,j)  =  0.0 
sc_byn2(i,j)  =  0.0 
sc_bxnh3(i,j)  =  0.0 
sc_bynh3(i, j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc^ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
130  continue 

do  140  1=1, nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 
140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if(sc_ax(i, j) .ge.O. )then 
sc_bxn2 ( i , j ) *dn2 ( i , j ) 
sc_bxh2 ( i , j ) =dh2 ( i , j ) 
sc_bxnh3(i, j )=dnh3(i, j ) 
else 

sc_bxn2 ( i , j ) =dn2 ( i+1 , j ) 
sc“bxh2(i, j)*dh2(i+l, j) 
s c~bxnh3 ( i , j ) =dnh3 ( i + 1 , j ) 
endi? 

150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if<sc_ay(i, j).ge.0. )then 
sc”byn2(i,j)  =  dn2(i,j) 
sc_byh2(i,j)  =  dh2(i,j) 
sc_bynh3(i, j)  =  dnh3(i,j) 
else  ~ 

sc_byn2(i,j)  =  dn2(i,j+l) 
sc“byh2(i,j)  =  dh2(i,j+l) 
sc”bynh3(i,j)  =  dnh3(i,j+l) 
endif~ 

160  continue 
c 
c 

do  170  i>l,nxl 
do  170  j=i,ny 

fh2(i,j)  «  sc_ax(i,j)  *  sc_bxh2(i,j) 

fn2(i,j)  =  sc~ax(i,j)  *  sc_bxn2(i,j) 

fnh3(i,j)  *  sc_ax(i,j)  *  sc_bxnh3(i, j ) 

170  continue  ” 

c 

do  180  i-l,nx 
c 

do  180  j»l,nyl 

gh2(l,j)  -  sc_ay(i,j)  *  sc_byh2(l,j) 

gn2(i,j)  =  sc_ay(i,j)  *  sc_byn2(i,j) 

gnh3(i,j)  *  sc_ay(i,j)  *  sc_bynh3U»j) 


180  continue 

do  200  i=2,nx 
do  200  j=2,ny 

dh2(i,j)  =  dh2(i,j)  -dt*(fh2(i , j )-fh2(i-l, j ) )/dx 
&  -dt*(gh2(i,j)-gh2(i,j-l))/dy 

&  -dt*dh2(i,j)*(cl*dopl(i,j)  +c3*dohpl(i,j^ 

&  +  cph3h2*dph3(i,j)  +  cph4h2*dph4(i , j )  +  cph5h2*dph5(i,j)) 

if  (abs(vxn(i,j)-vxoxy(i,j))  .ge.  10.2e5)  then  iVcrit  =  10.2  km/s 

if  (dciv(i,j).ge.civn21  .and.  dciv(i , j ) . le. civn22)  then 
cl0=(5.e-7*dciv(i,j)  -  0.17)*(vcivn2(i,i)/1.5)**5 
cl0=cl0*omegaN2 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i,  j)-fn2(i-l,-;  ))/dx 
&  'dt*(gn2(i,j)-gn2(i,j-l))/dy 

&  -dt*dn2(i,j)*(c5*dopl(i,j) 

&  +  cph3n2*dph3(i,j)  +  cph4n2*dph4(i , j )  +  cph5n2*dph5(i, i) 

&  +  cph6n2*dph6(i, j)  +  cph7n2*dph7(i, j )) 

^  -de(i,j)*(exp(clO*dt)-l.) 

end  if 


& 

& 

& 

& 

& 


if  (dciv(i,j).gt.civn22  .and.  dciv(i , j ) . le. civn23)  then 
clO=(7.8  -  1.2e-7*dciv(i, j))*(vcivn2(i- ,)/1.5)**5 
cl0»cl0*omegaN2 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i,j)-fn2(i-l, j))/dx 

-dt*(gn2(i,j)_gn2(i,j-l))/dy 

•dt*dn2(i,j)*(c3*dopl(i,j) 

+  cph3n2*dph3(i,j)  +  cph4n2*dph4(i, j)  +  cph5n2*dph5(i, 
+  cph6n2*dph6(i,j)  +  cph7n2*dph7(i, j)) 

~de(i,j)*(exp(clO*dt)-l.) 

end  if 


j) 


else 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i,j)^fn2(i-l, j))/dx 
&  -dt*(gn2(i,j)-gn2(i,j-l))/dy 

&  -dt*dn2(i,j)*(c5*dopl(i,j) 

&  +  cph3n2*dph3(i,j)  +  cph4n2*dph4(i, j )  +  cph5n2*dph5(i, i) 

&  +  cph6n2*dph6(i,j)  +  cph7n2*dph7(i, j)) 

end  if 

dnh3(i,j)  =  dnh3(i , j )-dt*(fnh3(i , j )-fnh3{i-l, j ) )/dx 
t  ^  -dt*(gnh3(i,j)-gnh3(i,j-l))/dy 

&  -dt*dnh3(i, j)*(c7*dopl(i,j) 

&  +  cphlnh3*dphl(i,j)  +  cph2na3*dph2(i, j)  +  cph3nh3*dph3(i,j) 

&  +  cph4nh3*dph4(i,j)*(80./170.)) 


200  continue 
return 
end 


subroutine  neutralbi 

parameter  (nx=90,ny=20) 

common/gr id/dx , dy , d  t 

common/ f luxH2 / f h2 ( nx , ny ) , gh2 ( nx , ny ) 

common  /  f  luxN2  /  f  n2  ( nx ,  ny ) ,  gn2  ( nx ,  ny  ) 

common/fluxH20/fh2o(nx,ny),gh2o(nx,ny) 

common/ f luxC02 / f  co2 ( nx , ny ) , gco2 ( nx , ny ) 
common/ fluxCO/  f  co ( nx ,  ny ) ,  geo ( nx ,  ny ) 
common / dens / dh2 ( nx , ny ) , dn  2 ( nx , ny ) , dnh3 ( nx , ny ) , 
$  dh2o(nx,ny),dco(nx,ny),dco2(nx,ny), 


$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opl(nx,ny),dnopl(nx,ny) , 

$  dnh3pl ( nx , ny ) , dn2pl ( nx , ny ) , de ( nx , ny ) , dh2pl ( nx , ny ) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opi(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny) ,dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloneu/vxn(nx,ny) , vyn(nx,ny) 
common/ velo i on/ vxoxy ( nx , ny ) , vy oxy { nx , ny ) 
common/ rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/rxratebi/cll,cl2,ci3,cl4,cl5,ci6,cl7,cl8,cl9 
common/ scrat_a/s  c_ax ( nx , ny ) , s  c_ay ( nx , ny ) 
common / s  cr a  t _b / sc_bx ( nx , ny ) , s  c_by ( nx , ny ) 
common/scratn/sc_bxh2(nx,ny) ,sc_byh2(nx,ny) , 

$  sc_bxn2(nx,ny) ,sc_byn2(nx,ny) , 

$  sc_bxnh3(nx,ny) , sc_bynh3(nx,ny) , 

$  sc_bxh2o(nx,ny) ,sc_byh2o(nx,ny) , 

$  sc_bxco(nx,ny) ,sc_byco(nx,ny) , 

$  sc_bxco2(nx,ny) ,sc_byco2(nx,ny) , 

$  sc_bxphl(nx,ny) ,sc_byphl(nx,ny) ,sc_bxph2(nx,ny) , 

$  sc_byph2 ( nx , ny ) , sc_bxph3 (nx , ny ) , sc_byph3 (nx , ny ) , 

$  sc_bxph4(nx,ny) , sc_byph4(nx,ny) ,sc_bxph5(nx,ny) , 

$  sc_byph5(nx,ny) , sc_bxph6(nx,ny) ,sc_byph6(nx,ny) , 

$  s  c_bxph7 ( nx , ny ) , s  c_by ph7 ( nx , ny ) 
common/ci vn21im/civn21 , civn22 , civn23 
common/civco21im/civco21 , civco22, civco23 
common/civ/dciv(nx,ny) , vcivn2(nx,ny) , vcivco2(nx,ny) 
common/ omega/ omeg2iN2 ,  omegaC02 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  =0.0 
sc_bxh2(i,j)  =  0.0 
sc_byh2(i,j)  =  0.0 
sc_bxn2(i,j)  =  0.0 
sc_byn2(i,j)  =  0.0 
sc_bxh2o(i, j)  =  0.0 
sc_byh2o(i , j )  =  0.0 
sc_bxco(i,j)  =  0.0 
sc_byco(i,j)  =  0.0 
sc_bxco2(i , j )  =  0.0 
sc_byco2(i,j)  =  0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if (sc_ax(i , j ) .ge.O. )then 


sc_bxn2(i, j )=dn2(i, 3 ) 
sc_bxh2(i,j)=dh2(i, j) 
sc_bxh2o(i , j )=dh2o(i ,3 ) 
sc_bxco( i , j )=dco( i , j ) 
sc_bxco2( i ,3 )=dco2(i ,3 ) 
else 

sc_bxn2(i , j )=dn2(i+l,3 ) 
sc_bxh2(i ,3 )=dh2(i+l,3 ) 
sc_bxh2o( i , 3 )=dh2o( i+1 ,3 ) 
sc_bxco(i,3 )=dco(i+l, j ) 
sc_bxco2(i , j )=dco2(i+l, j ) 
endif 
150  continue 
c 

do  160  i=:l,nx 
do  160  3=1, nyl 
if (sc_ay(i, j ) .ge.O. )then 
sc_byn2(i,3)  =  dn2(i,3) 
sc_byh2(i,3)  =  dh2(i,j) 
sc_byh2o(i ,3 )  =  dh2o(i,3) 
sc_byco(i,3)  =  dco(i,j) 
sc_byco2(i, j)  =  dco2(i,3) 
else 

sc_byn2(i,3)  =  dn2(i,j+l) 
sc_byh2(i,3)  =  dh2(i,j+l) 
sc_byh2o(i , j )  =  dh2o(i,3+l) 
sc_byco(i,3)  =  dco(i,3+l) 
sc_byco2(i,3)  =  dco2(i,3+l) 
endif 
160  continue 
c 
c 

do  170  i=l,nxl 
do  170  3=1, ny 

fh2(i,3)  =  sc_ax(i,3)  *  sc_bxh2(i,j) 

fn2(i,j)  =  sc_ax(i,3)  *  sc_bxn2(i,j) 

fh2o(i,3)  =  sc_ax(i,j)  *  sc  bxh2o(i,j) 
fco(i,3)  =  sc_ax(i,j)  *  sc_Bxco(i,j) 
fco2(i,3)  =  sc_ax(i,3)  *  sc_bxco2(i,3 ) 

170  continue 
c 

do  180  i=l,nx 
c 

do  180  3=1, nyl 

gh2(i,j)  =  sc_ay(i,3)  *  sc_byh2(i,j) 

gn2(i,3)  =  sc_ay(i,3)  *  sc_byn2(i,3) 

gh2o(i,3)  =  sc_ay(i,3)  *  sc  byh2o(i,3) 
gco(i,j)  =  sc_ay(i,3)  *  sc_Byco(i,3) 
gco2(i,j)  =  sc_ay(i,j)  *  sc_byco2(i,3 ) 

180  continue 
c 

do  200  i=2,nx 

do  200  j=2,ny 

dh2(i,j)  =  dh2(i,3)  -dt*(fh2(i,j)-fh2(i-l,3))/dx 
&  -dt*(gh2(i,j)-gh2(i,j-l))/dy 

&  -dt*dh2(i ,3 )*(cl*dopl(i, j )  +c3*dohpl(i,3)) 


if  (abs(vxn(i,3)-vxoxy(i,3)).ge.l0.2e5)  then 
if  (dciv(i,3).<je.civn21  .and.  dciv{i,j).lc.civn22)  then 


Q^  Q^  Q**  tt^  Q^  Q^ 


clO=(5.e-7*dciv(i, j)  -  0.17)*(vcivn2(i, j)/1.5)**5 
clO=clO*omegaN2 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i , j )-fn2(i-l, j ) )/dx 

-dt*(gn2(i, j)-gn2(i, j-l))/dy 
-d  t *dn2 ( i , j )*c5*dopl ( i , j ) 
-de(i,j)*(exp(clO*dt)-l. ) 

end  if 

if  (dciv(i, j).gt.civn22  .and.  dciv(i , j ) . le. civn23)  then 
clO=(7.8  -  1.2e-7*dciv(i, j))*(vcivn2(i,j)/1.5)**5 
clO=clO*omegaN2 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i,j)-fn2(i-l,j))/dx 

-dt*(gn2(i,j)-gn2(i, j-l))/dy 
-d  t *dn2 ( i , j ) *c5*dopl ( i , j ) 
-de(i,j)*(exp(clO*dt)-i. ) 

end  if 
else 

dn2(i,j)  =  dn2(i,j)  -dt*(fn2(i,j)-fn2(i-l, j))/dx 

-dt*(gn2(i,j)-gn2(i,j-l))/dy 
-dt*dn2(i, j)*c5*dopl(i, j) 

end  if 

dh2o(i,j)  =  dh2o(i, j)-dt*(fh2o(i,j)-fh2o(i-l,j))/dx 

-dt*(gh2o(i,j)-gh2o(i, j-l))/dy 
-dt*(dh2o( i , j )*cll*dopl( i , j ) ) 

-dt*(dh2o(i, j )*cl8*dh2opl(i, j)) 

dco(i,j)  =  dco(i, j)-dt*(fco(i,j)-fco(i-l,j))/dx 

-dt*(gco(i,j)-gco(i,j-lb/dy 
+dt*(dco2(i, j )*cl2*dopl(i, j )) 

if  (abs(vxn(i, j)-vxoxy(i, j))  .ge.  7.7e5)  then 

if  (dciv(i, j).ge.civco21  .and.  dciv(i, j).le.civco22)  then 
cl3=(0.25  +  4.5e-7*dciv(i, j))*(vcivco2(i, j)/1.5)**5 
cl3=cl3*omegaC02 

dco2(i,j)  =  dco2(i, j )-dt*(fco2(i,j)-fco2(i-l, j ))/dx 

-dt*(gco2(i,j)-gco2(i,3-lh/dy 
-dt*(dco2(i , j )*cl2*dopl(i , j ) ) 

-dt*(dco2(i, j )*cl4*dopl(i, j )) 

-de(i, j )*(exp(cl3*dt)-l. ) 

end  if 

if  (dciv(i,j).gt.civco22  .and.  dciv(i,j).le.civco23)  then 
cl3=(17.  -  2.23e-7*dciv(i,j))*(vcivco2(i,j)/1.5)**5 
cl3=cl3*omegaC02 

dco2(i,j)  =  dco2(i, j)-dt*(fco2(i,j)-fco2(i-l,j))/dx 

-dt*(gco2(i,j)-gco2(i,j-l))/dy 
-dt*(dco2(i , j )*cl2*dopl(i, j )) 

-dt*(dco2(i, j )*cl4*dopl(i, j)) 

-de(i, j )*(exp(cl3*dt)-l. ) 

end  if 
else 

dco2(i,3)  *  dco2(i,3)-dt*(fco2(i,j)-fco2(i-l,3))/dx 

-dt*(gco2(i,3)-gco2(i,j-l))/dy 
-dt*(dco2(i ,3 )*cl2*dopl(i,3) ) 

-dt*(dco2(i, j)*cl4*dopl(i, j)) 


endif 


200  continue 
return 
end 
c 

subroutine  oplus 
c 

parameter  (nx=90,ny=20) 
common/gr id/dx , dy , d  t 

common/ f luxOpl/ f opl ( nx , ny ) , gopl( nx , ny ) 
common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) , 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opl(nx,ny) ,dnopl(nx,ny) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opi(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny) ,dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
c  ommon /veloion/vxoxy(nx,ny), vy oxy ( nx , ny ) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common  /  s  c  ra  t_a/  s  c_ax  ( nx ,  ny ) ,  s  c_ay  (  nx ,  ny  ) 
common/scrat_b/sc_bx(nx,ny) ,sc_by(nx,ny) 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=i,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  =  0.0 
sc_bx(i,j)  =  0.0 
sc_by(i,j)  »  0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

c  vrx  =sc_ax(i,j)=0.5*(vx(i+l,j)+vx(i,j)) 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxoxy(i+l, j)  +  vxoxy(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyoxy(iJ+l)  +  vyoxy(i,j)) 

140  continue 
c 

do  150  i*l,nxl 
do  150  j*l,ny 

if (sc_ax(i, j).ge.0. )then 
sc_bx(i,j)=dopl(i,j) 
else 

sc  bx(i, j)*dopl(i+l,j) 
endil 

150  continue 
c 

do  160  i*l,nx 
do  160  j=l,nyl 
i f ( sc_ay ( i , j ) . ge . 0 . ) then 
sc_by(i,j)  -  dopl(i,j) 
else 


sc_by(i,j)  =  (iopl(i,j+l) 

end  if 

160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=l,ny 

fopl(i,j)  =  sc_ax(i,j)  *  sc_bx(i,j) 
170  continue 


do  180  i=l,nx 
do  180  j=l,nyl 

gopl(i,j)  =  sc_ay(i,j)  *  sc  by(i,j) 
180  continue 


c 


200 


do  200  i=2,nx 
do  200  j=2,ny 

dopl(i,j)  =  dopl(i,j)  -dt*(fopl(i,j)-fopl(i-l,j))/dx 
f  j  .  -<it*(gopl(i,j)-gopl(i,j-l))/dy 

&  -dt*dopl(i,j)*(dh2(i,j)*cl  +  dn2(i,j)*c5  +  dnh3(i , j )*c7) 
continue 
return 
end 


c 

subroutine  oplusbi 
c 


parameter  (nx=90,ny=20) 

common/grid/dx,dy,dt 

common/ f luxOpl/ f opl ( nx , ny ) , gopl ( nx , ny ) 

common/dens /dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 (nx , ny ) , 

S  an2o(nx,ny),dco(nx,ny),dco2(nx,ny), 

$  dopl(nx,ny),dohpl(nx,ny),dh2opl(nx,ny),dnopl(nx,ny), 

S  dnh3pl(nx,ny) »dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 

$  ao2pl(nx,ny),dcopl(nx,ny),dco2pl(nx,ny),dh3opl(nx,ny), 
$  dphl(nX)ny) >dph2(nx,ny) ,dph3(nx,ny) ,dph4(nx,ny) , 

$  dph5(nx,ny),dph6(nx,ny),dph7(nx,ny) 
commou/veloion/vxoxy(nx,ny),vyoxy(nx,ny) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/rxratebi/cll , cl2 , cl3 , cl4 , cl5 , cl6 , cl7 , cl8 , cl9 
common / s era t  a/ s c_ax ( nx , ny ), s cay ( nx , ny ) 

common/ scrat_b/sc_bx(nx,ny),s  c_by ( nx , ny ) 


-  initialize  scratch  arrays 

do  100  i=l,nx 
do  100  jsl,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  *  0.0 

sc~bx(i»j)  =  0.0 

sc_by(i,j)  =  0.0 
100  continue 

nxl*nx-l 

nyl«ny-l 

do  130  i«l,nxl 
do  130  j»l,ny 

sc_ax(i,j)  s  0.5  *  (vxoxy(i+l,j)  +  vxoxy(i,i)) 
130  continue 

do  140  isl,nx 


do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyoxy(i, j+1)  +  vyoxy(i,j)) 

140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if (sc_ax(i , j ) .ge.O. )then 
sc_bx(i,j)=dopl(i, j) 
else 

sc_bx(i, j)*dopl(i+l,j) 
endif 
150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if (sc_ay(i , j ) .ge.O. )then 
sc_by(i,j)  =  dopl(i,j) 
else 

sc_by(i,j)  =  dopl(i,j+l) 
endif 
160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=l,ny 

fopl(i,j)  =  sc_ax(i,j)  *  sc_bx(i,j) 

170  continue  ~ 

c 

do  180  islynx 
do  180  j=l,nyl 

gopl(i,j)  *  sc_ay(i,j)  *  sc_by(i,j) 

180  continue 
c 

do  200  i=2,nx 
do  200  j=2,ny 

dopl(i,j)  =  dopl(i,j)  -dt*(fopl(i,j)-fopl(i-l,j))/dx 
&  -dt*(gopl(i,j)-gopl(i,j-l))/dy 

&  -dt*dopl(i,j)*(dh2(i, j)*cl  +  dn2(i,j)*c5  +  dh2o(i, j)*cll 
&  +  dco2(i , j )*cl2  +  dco2(i , j )*cl4  +  dco(i, j )*cl7) 

200  continue 
return 
end 
c 


subroutine  ohplus 
c 

parameter  (nx*90,ny=20) 
common/gr id/dx , dy , dt 

common/ f luxOHpl/ f ohpl ( nx , ny ) , gohpl ( nx , ny ) 
common/dens /dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny),dohpl(nx,ny),dh2opi(nx,ny),dnopl(nx,ny), 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny) , 

$  do2pl (nx , ny ) , dcopl ( nx , ny ) , dco2pl ( nx , ny ) , dhSopl ( nx , ny ) , 
$  dphl(nx,ny),dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny), 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/veloneu/ vxn ( nx , ny ) , vyn ( nx , ny ) 
common/rxratc/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9, clO 
common/rxratebi/cll,cl2,cl3,cl4,cl5,cl6,cl7,cl8,cl9 


common/scrat_a/sc_ax(nx,ny) ,sc  ay(nx,ny) 
common/scrat_b/sc_bx(nx,ny),sc_by(nx,ny) 


c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  =0.0 

sc_bx(i,j)  =  0.0 

sc_by(i,j)  =  0.0 
loo  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 
140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if(sc_ax(i, j),ge.0. )then 
sc~bx(i,  j)«dohplU,j) 
else 

sc  bx(i , j )«dohpl(i+l, j ) 
endil 

150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if(sc_ay(i,j).ge.0. )then 
sc_by(i,j)  =  dohpl(i,j) 
else 

sc_by(i,j)  =  dohpl(i,j+l) 
end  if 

160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=i,ny 

fohpKiJ)  =  sc_ax(i,j)  *  sc_bx(i,j) 

170  continue 


c 

do  180  ix;l,nx 
do  180  j»l,nyl 

gohpl(i,j)  =  sc_ay(i,j)  *  sc_by(i,j) 
180  continue 


do  200  i=2,nx 
do  200  j«2,ny 

dohpl(i,j)  =  dohpl<i,j)  -dt*(fDhpl(i,j)-fohpl(i-l,j))/dx 
&  -dt*(gohpl(i,j)-gohpl(i,j-l))/dy 

&  ^dt*(dopl(i,j)*cl*dh2(l,j)  -  dohpKi, j)*c2*de(i, i) 

6  -  dohpl(i, j)*c3*dh2(i,j)) 


200  continue 
return 
end 
c 


subroutine  hZoplus 

parameter  (nx=90,ny=20) 
common/grid/dx , dy , d t 

common/ fluxH20pl/ fh2opl ( nx , ny ) , gh2opl (nx , ny ) 
common/dens /dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 ( nx , ny ) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) , 

$  dopl ( nx , ny ) , dohpl ( nx , ny ) , dh2opl ( nx , ny ) , dnopl ( nx , ny ) , 

$  dnh3pl(nx,ny),dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny), 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opl(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/veloneu/vxn(nx,ny) ,vyn(nx,ny) 
common/veloion/vxoxy(nx,ny),vyoxy(nx,ny) 
common/rxra te/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common /scrata/s  c_ax ( nx , ny ) , s  cay ( nx , ny ) 

c ommon / s c r a t _b/ s c_bx ( nx , ny ) , s c_by ( nx , ny ) 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=l,ny 
sc  ax(i,j)  *  0.0 
sc"ay(i,j)  =  0.0 
sc_bx(i,j)  =  0.0 
sc“by(i,j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc^ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  ial,nxl 
do  150  j=l,ny 
if(sc_av(:,j).ge.0.)then 
sc_bx(i,j).dh2opl(i,j) 
else 

sc  bx(i, j)*dh2opl(i+l, j) 
endil 

150  continue 
c 

do  160  i=l,nx 
do  160  j»l,nyl 
if(sc_ay(i,3).ge.0. )then 


160 

c 

c 


170 

c 


180 


sc_by(i,j)  =  dh2opl(i,j) 

else 

sc_by(i,j)  =  dh2opl(i , j+1) 
end if” 
continue 


do  170  i=l,nxl 
do  170  j=l,ny 

fh2opl(i,j)  =  sc_ax(i,j)  *  sc_bx(i,j) 
continue 

do  180  i=l,nx 
do  180  j=l,nyl 

gh2opl(i,j)  =  sc_ay(i,j)  *  sc_by(i,j) 
continue 


do  200  i=2,nx 
do  200  j=2,ny 

dh2opl(i,j)  =  dh2opl(i,j)  -dt*(fh2opl(i,j)-fh2opl(i-l,j))/dx 
&  -dt*(gh2opl(i,j)-gh2opl(i,j-l))/dy 

&  +dt*(dohpl(i, j)*c3*dh2(i, j)  -  dh2opl(i, j)*c4*de(i,j)) 

200  continue 
return 
end 


c 


subroutine  h2oplusbi 
c 


parameter  (nx=90,ny*20) 
common/grid/dx , dy , d t 

common/ fluxH20pl/fh2opl(nx,ny) ,gh2opl(nx,ny) 
common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opi(nx,ny),dnopl(nx,ny) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny),dh3opl(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
c  ommon  /  ve  loneu  /  vxn  (  nx ,  ny ) ,  vy n  (  nx ,  ny  ) 
common/ rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/rxratebi/cll,cl2,ci3,cl4,cl5,cl6,cl7,cl8,cl9 
c ommo n / s c r a t _a / s c_ax ( nx , ny ) , s c_ay ( nx , ny ) 
common/ s  c  ra  t_b/ s  c_bx ( nx , ny ) , s  c_by ( nx , ny ) 

c -  initialize  scratch  arrays 

c 

do  100  i«l,nx 
do  100  j>il,ny 
sc_ax(i,j)  *  0.0 
sc_ay(i,j)  *  0.0 
sc_bx(i,j)  «  0.0 
sc_by(i,j)  *  0.0 
100  continue 
c 

nxl*nx-l 

nyl»ny-l 

c 

do  130  i-l,nxl 


do  130  j=l,ny 

sc_ax(i,i)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 
140  continue 


c 

do  150  i=l,nxl 
do  150  j=l,ny 

if(sc_axa,j).ge.0.)then 
.sc_bx(  i ,  j  )=dh2opl(  i ,  j  ) 
else 

sc  bx(i,j)*dh2opl(i+l,j) 
endil 
150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if(sc_ay(i, j).ge.0. )then 
sc_by(i,j)  =  dh2opl(i,j) 
else 

sc_by(i,j)  =  dh2opl(i,j+l) 
end  if 
160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=i,ny 

fh2opl(i,j)  »  sc  ax(i,3)  *  sc  bx(l,j) 
170  continue  ~  ~ 


do  180  isl,nx 
do  180  ja:l,nyl 

gh2opl(i,j)  «  sc_ay(i,j)  *  sc  by(i,j) 
180  continue 


200 


do  200  i32,nx 
do  200  js2,ny 

dh2opl(i,j)  =  dh2opl(i,j)  -dt*(fh2opl(i, j)-fh2opl(i-l,j))/dx 
^  ,  -dt*(gh2opl(i,j)-gh2opl(i,j-l))/dy 

+dt*(dohpl(i,j)*c3*dh2(i,j)  -  dh2opl(i, j)*c4*de(i, j) 

+  dh2o(i,j)*cll*dopl(i,j)  -  dh2opl(i,j)*cl8*dh2o(i,i)) 
continue 
return 
end 


c 


subroutine  noplus 

parameter  (nx»90,ny»20) 
common/gr id/dx , dy , d  t 

common/ £ luxNOpl / f nopl ( nx , ny ) , gnopl ( nx , ny ) 
common/dens/dh2 (nx , ny ) , dn2 ( nx , ny ) , dnh3 (nx , ny ) , 

$  dh2o(nx,ny),dco(nx,ny),dco2(nx,ny), 

$  dopl( nx , ny ) , dohpl ( nx , ny ) , dh2opl(nx , ny ) , dnopl (nx , ny ) , 

$  dnh3pl(nx,ny),dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny), 

$  do2pl ( nx , ny ) , dcopl ( nx , ny ) , dco2pl ( nx , ny ) , dh3opl ( nx , ny ) , 
$  dphl(nx,ny),dph2(nx,ny),dph3(nx,ny),dph4(nx,ny), 


$  dph5(nx,ny),dph6(nx,ny),dph7(nx,ny) 
common/ ve loneu/vxn ( nx , ny ) , vyn ( nx , ny ) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 

common/rxratebi/cll,cl2,cl3,cl4,cl5,cl6,cl7,cl8,cl9 
common / s c r a t_a / s c_ax ( nx , ny ), s c_ay ( nx , ny ) 

common/ s  c  r a  t_b/ s  c_bx ( nx , ny ) , sc_by ( nx , ny ) 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =  0.0 
sc_ay(i,j)  =  0.0 

sc_bx(i,j)  =  0.0 

sc_by(i,j)  =  0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

=  0*5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

1/A  =  0-5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  i=l,nxl 
do  150  j»l,ny 
if(sc_ax(i,j).ge.0.)then 
sc_bx(i, j)=dnopl(i, j ) 
else 

sc  bx<i,j)»dnopl(i+l,j) 
endi? 

150  continue 
c 

do  160  i>l,nx 
do  160  j=l,nyl 
if(sc_ay(i,j).ge.0.)then 
sc_by(i,j)  *  dnopl(i,j) 
else 

sc_by(i,j)  *  dnopl(i,j+l) 
end  if 

160  continue 
c 
c 

do  170  i*l,nxl 
do  170  j=i,ny 

i-,A  *  sc  ax(i,j)  *  sc  bx(i,j) 

170  continue  ”  ~ 

c 

do  180  i«l,nx 
do  180  j»l,nyl 

10A  =  sc_ay(i,j)  *  sc  by(i,j) 

180  continue  ~ 

c 

do  200  i*2,nx 
do  200  j«2,ny 


dnopl(i,j)  =  dnopKiJ)  -dt*(fnopl(i,j)-fnopl(i-l, j))/dx 
£,  -dt*(gnopl(i,j)-gnopl(i,j-l))/dy 

&  +dt*(dopl(i,j)*c5*dn2(i,j)  -  dnopl(i,j)*c6*de(i,j)) 
200  continue 
return 
end 

subroutine  nh3plus 

c 

parameter  (nx=90,ny=20) 
common/grid/dx , dy , d t 

common/ £luxNH3pl/ fnh3pl ( nx , ny ) , gnh3pl ( nx , ny ) 
common/dens/dh2(nx,ny)  ,<in2(nx,ny)  ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny) , 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opl(nx,ny) ,dnopl(nx,ny) , 

$  dnh3pl(nx,ny),dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny), 

$  do2pl(nx,ny),dcopl(nx,ny),dco2pl(nx,ny),dh3opl(nx,ny), 

$  dphl(nx,ny),dph2(nx,ny),dph3(nx,ny),dph4(nx,ny), 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/veloneu/vxn(nx,ny),vyn(nx,ny) 
common/ veloph/vxph ( nx , ny ) , vyph ( nx , ny ) 
common/ rxrate/ cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/phrateH2/cph3h2 , cph4h2 , cph5h2 
common/phrateN2/cph3n2 , cph4n2 , cph5n2 , cph6n2 , cph7n2 
common/phrateNH3/cphlnh3 , cph2nh3 , cph3nh3 , cph4nh3 
common/scrat  a/sc_ax(nx,ny) ,sc_ay(nx,ny) 
common/scrat”b/sc”bx(nx,ny),sc'by(nx,ny) 

c -  initialize  scratch  arrays 

c 

do  100  i*l,nx 
do  100  j«l,ny 
sc_ax(i,j)  =0.0 
sc“ay(i,j)  =  0.0 
sc”bx(i,j)  =0.0 
sc”by(i,j)  =0.0 
100  continue 
c 

nxl*nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  3=1, nyl 

sc  ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  i*l,nxl 
do  150  j«l,ny 

if(sc  ax(i,j).ge.0. )then 
sc'bx ( i , j ) =dnh3pl ( i , j ) 
else  ~ 

sc  bx(i,j)«dnh3pl(i+l,j) 
endil 

150  continue 
c 

do  160  i=l,nx 


do  160  j=l,nyl 
i f ( sc_ay ( i , j ) . ge . 0 . ) then 
sc_by(i,j)  =  dnh3pl(i,j) 
else 

sc_by(i,j)  =  dnh3pl(i, j+1) 
end  if 

160  continue 
c 


do  170  i=l,nxl 
do  170  j=i,ny 

fnh3pl(i,j)  =  sc_ax(i,j)  *  sc_bx(i,j) 

170  continue 
c 

do  180  i=l,nx 
do  180  j=l,nyl 

gnh3pl(i,j)  =  sc_ay(i,j)  *  sc_by(i,j) 

180  continue  ~ 

c 

do  200  i=2,nx 
do  200  j=2,ny 

dnh3pl(i,j)  =  dnh3pl(i,j)  -dt*(fnh3pl(i, j )-fnh3pl(i-l, j ))/dx 
&  -dt*(gnh3pl(i,j)-gnh3pl(i,j-l))/dy 

&  +dt*(dopl(i,j)*c7*dnh3(i, j)  -  dnh3pl(i, j)*c8*de(i,j)) 

&  +dt*dnh3(i,j)*(cphlnh3*dphi(i,j)  +  cph2nh3*dph2(i , j ) 

&  +cph3nh3*dph3(i, j )  +  cph4nh3*dph4(i,j)*(80./170. )) 

200  continue 
return 
end 
c 

subroutine  n2plus 
c 

parameter  (nx=90,ny=20) 
common/grid/dx,dy ,dt 

comfflon/fluxN2pi/£n2pl(nx,ny),gn2pl(nx,ny) 
common/dens /dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 ( nx , ny ) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opl(nx,ny),dnopl(nx,ny) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny),dco2pl(nx,ny),dh3opl(nx,ny) , 

$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny), 

$  dph5(nx,ny),dph6(nx,ny),dph7(nx,ny) 
common  /  ve  loneu  /  vxn  (  nx ,  ny ),  vyn  (  nx ,  ny  ) 
common/ veloph/vxph ( nx , ny ) , vyph (nx , ny ) 
common /veloion/vxoxy(nx,ny), vy oxy ( nx , ny ) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , ciO 
common/rxratebi/cll,cl2,ci3,cl4,cl5,ci6,cl7,cl8,cl9 
common/phrateN2/ cph3n2 , cph4n2 , cph5n2 , cph6n2 , cph7n2 
common/scrat_a/sc_ax(nx,ny),sc_ay(nx,ny) 
common/scrat  b/sc  bx(nx,ny) ,sc_by(nx,ny) 
common/civn2Tim/cIvn21 , ci vn22 , civn23 
common/civ/dciv(nx,ny) , vcivn2(nx,ny) ,vcivco2(nx,ny) 
cofflmon/omega/omegaN2 , omegaC02 

common/elrestl/elfluence, in2restrict,ico2restrict 
coinmon/elrest2/dn2plcheck(nx,ny) ,dco2plcheck(nx,ny) 

c -  initialize  scratch  arrays 

c 


do  100  i=l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =0.0 

sc_ay(i,j)  =0.0 

sc_bx(i,j)  =  0.0 
sc_by(i,j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc  ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
130  ''op.tinue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  wn(i,j)) 
140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if(sc_ax(i,j).ge.0.)then 

sc_bx(i,j)=dn2pl(i,j; 

else 

sc  bx(i, j)=dn2pl(i+l, j) 
endiT 

150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if(s''_ay(i,j).gc.0.  )then 
sc_by(i,j)  =  dn2pl(i,j) 
else 

sc_by(i,j)  =  dn2pl(i,j+l) 
end  if 

160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=l,ny 

fn2pni,j)  =  sc_ax(i,j)  *  sc  bx(i,j) 

170  continue 
c 

do  180  i=l,nx 
do  180  j=l,nyl 

gn2pl(i,j)  =  sc_ay(i,j)  *  sc_by(i,j) 

180  continue 
c 

do  200  i=2,nx 
do  200  j=2,ny 

if  (abs(vxn(i,j)-vxoxy(i,j))  .ge.  10.2e5)  then 


& 

& 

& 


if  (dciv(i,j).ge.civn21  .and.  dciv(i,j).le.civn22)  then 
clO=(5.e-7*dciv(i,j)  -  0. 17)*<vcivn2(i, j )/1.5)**5 
cl0=cl0*omegaN2 

dn2pl(i,j)  .  dnrpl(i,j)  -dt*(fn2pl(i,j)-fn2pl(i-l, j))/dx 

-dt*(gn2pl(i,j)-gn2pl(i,j-l))  'dy 
+dt*dn2(i,j)*(cph3n2*dph3(i,j)  +  cph4n2*dph4(i , j ) 
t-  cph5n2*dph5(i,j)  +  cph6n2*dph6(i,j)  +  cph7n2*dph7( ij)) 


«?»<?>  O' 


&  +de(i,j)*(exp(clO*dt)-l.) 

deltan2pl  =  dn2pl(i,j)  -  dn2plcheck,(i, j ) 
if  (deltan2pl.gt.elfluence)  then 
dn2pl(i,j)  =  dn2plcheck(i , 3 )  +  elfluence 
in2restrict  =  in2restrict  +  1 
end  if 

dn2plcheck(i, j )=dn2pl(i, j ) 
end  if 

if  (dciv(i , j ) .gt . civn22  .and.  dciv(i,3).le.civn23)  then 
clO=(7.8  -  1.2e-7*dciv(i,j))*(vcivn2(i,j)/1.5)**5 
clO=clO*omegaN2 

dn2pl(i,j)  =  dn2pl(i,j)  -dt*(fn2pl(i , j )-fn2pl(i-l, j ) )/dx 

-dt*(gn2pl(i,j)-gn2pl(i,j-l))/dy 
+dt*dn2(i , j )*(cph3n2*dph3(i,j )  +  cph4n2*dph4(i ,3 ) 

+  cph5n2*dph5(i,3 )  +  cph6n2*dph6(i, j )  +  cph7n2*dph7(i,3)) 
+de(i,3 )*(exp(clO*dt)-l. ) 
deltanipl  =  dn2pl(i,j)  -  dn2plcheck(i, j) 
if  (deltan2pl.gt. elfluence)  then 
dn2pl(i,j)  =  dn2plcheck(i, j )  +  elfluence 
in2restrict  =  in2restrict  +  1 
end  if 

dn2plcheck( i , j )=dn2pl( i , j ) 
end  if 

else 

dn2pl(i,j)  =  dn2pl(i,3)  -dt*(fn2pl(i,3 )-fn2pl(i-lo ))/dx 
&  -dt*(gn2pl(i,3)-gn2pl(i,3-l))/dy 

&  +dt*dn2(i, j)*(cph3n2*dph3(i,3)  +  cph4n2*dph4(i, j) 

&  +  cph5n2*dph5(i, j )  +  cph6n2*dph6(i, j )  +  cph7n2*dph7(i, j)) 

end  if 

200  continue 
return 
end 


subroutine  h2plus 

parameter  (nx=90,ny=20) 
common/gr id/dx , dy , d  t 

common/fluxH2pi/fh2pl(nx,ny) ,gh2pl(nx,ny) 
common/dens /dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 ( nx , ny ) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opl(nx,ny) ,dnopl(nx,ny) , 

$  dnh3pl ( nx , ny ) , dn2  pi ( nx , ny ) , de ( nx , ny ) , dh2pl ( nx , ny ) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny),dh3opl(nx,ny) , 
$  dphl(nx,ny),dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/veloneu/vxn(nx,ny) ,vyn(nx,ny) 
common/veloph/vxph(nx,ny) ,vyph(nx,ny) 
common/veloion/vxoxy(nx,ny) ,vyoxy(nx,ny) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/phrateH2/cph3h2 , cph4h2 , cphSh2 
common/scrat_a/sc_ax(nx,ny),sc_ay(nx,ny) 
common / s c r a t_b/ s c_bx ( nx , ny ), s c_by < nx , ny ) 

initialize  scratch  arrays 


do  100  i=l,nx 
do  100  j=l,ny 

sc_ax(i,j)  =  0.0 

sc_ay(i,j)  =  0.0 
sc_bx(i,j)  =  0.0 
sc_by(i,j)  =  0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130.  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 
140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 
if(sc_ax(i, j).ge.0. )then 
sc_bx(i,j)=dh2pl(i,j) 
else 

sc  bx(i,j)=dh2pl(i+l,j) 
endif 

150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if(sc  ay(i,j).ge.0.)then 
sc_by(i,j)  =  dh2pl(i,j) 
else 

sc_by(i,j)  =  dh2pl(i,j+l) 
end  if 

160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j=l,ny 

fh2pl(i,j)  =  sc_ax(i,j)  *  sc  bx(i,j) 

170  continue  ” 

c 

do  180  i*l,nx 
do  180  j=l,nyl 

gh2pl(i,j)  =  sc  ay(i,j)  *  sc  by(i,j) 

180  continue 
c 

do  200  i-2,nx 
do  200  j»2,ny 


& 

& 

& 


dh2pl(i,j)  »  dh2pl(i,j)  -dt*(fh2pl(i,j)-fh2pl(i-l,j))/dx 
.  -dt*(8h2pl(i,j)-gh2pl(i,j-l))/dy 

+dt*dh2(i,j)*(cph3h2*dph3(i,j)  +  cph4h2*dph4(i , j ) 

+  cph5h2*dph5(i,j)) 


200  continue 
return 
end 


c 


subroutine  co2plusbi 


parameter  (nx=90,ny=20) 
common/gri<l/dx,dy,dt 

common/fluxC02pl/fco2pl(nx,ny),gco2pl(nx,ny) 
common/ dens / dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 ( nx , ny ) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opi(nx,ny),dnopl(nx,ny) , 

$  dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny),dh2pl(nx,ny) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opi(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny) ,dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
c ommon / ve 1 on eu / vxn ( nx , ny ) , vyn ( nx , ny ) 
c ommon / ve 1 o i on / vxoxy ( nx , ny ) , vy oxy ( nx , ny ) 
common/rxrate/ cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/rxratebi/cll , cl2 , ci3 , cl4 , cl5, clS, cl7 , cl8 , cl9 
common/phrateC02/cph2co2 , cph3co2 , cph4co2 , cph5co2 
common/scrat_a/sc_ax(nx,ny) ,sc_ay(nx,ny) 
common/ s  c  ra  t_b / sc_bx ( nx , ny ) , s  c_by ( nx , ny ) 
comraon/ci vco21im/civco2i , civco22 , civco23 
common/civ/dciv(nx,ny) ,vcivn2(nx,ny) ,vcivco2(nx,ny) 
common/omega/omegaN2 , omegaC02 

common/elrestl/elfluence,in2restrict,ico2restrict 
common/elres  1 2  /dn2plcheck(  nx ,  ny  ) ,  dco2plcheck.(  nx ,  ny ) 

c -  initialize  scratch  arrays 

c 

do  100  isl,nx 
do  100  j=l,ny 
sc_ax(i,j)  =  0.0 
sc“ay(i,j)  =  0.0 
sc_bx(i,j)  =  0.0 
sc_by(i,j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  isl,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  i=l,nxl 
do  150  j«l,ny 

i f ( sc_ax( i , j ) . ge . 0 . ) then 

sc_bx(i, j)=dco2pl(i, j) 

else 

sc  bx(i,j)»dco2pl(i+l,j) 
endi? 

150  continue 
c 

do  160  i«l,nx 
do  160  j=l,nyl 


if(sc_ay(i,j).ge.O. )then 
sc_by(i,j)  =  dco2pl(i,j) 
else 

sc_by(i,j)  =  dco2pl(i, j+1) 
end  if 
160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j*i,ny 

fco2pl(i,j)  =  sc_ax(i,j)  *  sc_bx(i,j) 

170  continue 
c 

do  180  i=l,nx 
do  180  j=l,nyl 

gco2pl(i,j)  =  sc_ay(i,j)  *  sc_by(i,j) 

180  continue 
c 

do  200  i=2,nx 
do  200  j=2,ny 

if  (abs(vxn<i,j)-vxoxy(i,j))  .ge.  7.7e5)  then 

if  (dciv(i,j).ge.civco21  .and.  dciv(i, j).le.civco22)  then 
cl3=(0.25  +  4.5e-7*dciv(i,j))*(vcivco2(i, j)/1.5)**5 
cl3scl3*omegaC02 

dco2pl(i,j)  =  dco2pl(i,j)  -dt*(fco2pl(i,j)-fco2pl(i-l,j))/dx 
&  -dt*(gco2pl(i,j)-gco2pl(i, j-l))/dy 

&  +dt*(dopl(i, j)*cl4*dco2(i, j)  -  dco2pl(i, j)*cl5*de(i,j)) 

5  +de(i, j )*(exp(cl3*dt)-l. ) 

deltaco2pl  *  dco2pl(i,j)  -  dco2plcheck(i,j) 
if  (deltaco2pl.gt.elfiuence)  then 
dco2pl(i,j)  =  dco2plcheck(i,j)  elfluence 
ico2restrict  =  ico2restrict  +  1 
end  if 

dco2plcheck( i , j )»dco2pl( i , j ) 
end  if 

if  (dciv(i,j).gt.civco22  .and.  dciv(i, j).le.civco23)  then 
cl3*(17.  -  2.23e-7*dciv(i, j))*(vcivco2(i,j)/1.5)**5 
cl3=cl3*omegaC02 

dco2pl(i,j)  =  dco2pl(i,j)  -dt*(fco2pl(i,j)-fco2pl(i-l,j))/dx 

6  -dt*{gco2pl(i,  j)-gco2pl(i,  j-l))''<ly 

5  +dt*(dopl(i,j)*cl4*dco2(i,j)  -  dco2pl(i,j)*cl5*de(i,j)) 

6  +de(i, j)*(exp(cl3*dt)-l. ) 

c  Checking  ambient  electron  flux  available  to  neutralize  plasma 

deltaco2pl  *  dco2pl(i,j)  -  dco2plcheck(i, j) 
if  (deltaco2pl.gt. elfluence)  then 
dco2pl(i,j)  =  dco2plcheck(i, j)  +  elfluence 
ico2restrict  »  ico2restrict  +  1 
endif 

dco2plcheck( i , j )«dco2pl{ i , j ) 
endif 

else 

dco2pl(i,j)  .  dco2pl(i,j)  -dt*<fco2pl(i, j)-fco2pl(i-l, j))/dx 
&  -dt*(gco2pl(i,j)-gco2pl(i,j-l))/dy 

&  +dt*(dopl(i, j )*cl4*dco2(i, j )  -  dco2pi(i,j)*cl5*de(i,j)) 
endif 


200  continue 
return 
end 


subroutine  o2plusbi 
c 

parameter  (nx=90,ny=20) 
common/grid/dx , dy , dt 

common/flux02pi/fo2pl(nx,ny),go2pl(nx,ny) 
common/dens/dh2(nx,ny) ,dn2(nx,ny) ,dnh3(nx,ny) , 

$  dh2o(nx,ny) ,dco(nx,ny) ,dco2(nx,ny), 

$  dopl(nx,ny) ,dohpl(nx,ny) ,dh2opi(nx,ny),dnopl(nx,ny) , 

$  dnh3pl ( nx , ny ) , dn2pl ( nx , ny ) , de ( nx , ny ) , dh2pl ( nx , ny ) , 

$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny),dh3opi(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloneu/vxn( nx , ny ) , vyn ( nx , ny ) 
common/rxrate/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 
common/rxra tebi/cll , cl2 , ci3 , cl4 , cl5 , eld , cl7 , cl8 , cl9 
common/scrat_a/sc_ax(nx,ny) ,sc_ay(nx,ny) 
common / s c r a t”b / s c_bx ( nx , ny ), s c_by ( nx , ny ) 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j»l,ny 
sc_ax(i,j)  =0.0 
sc_ay(i,j)  «  0.0 
sc_bx(i,j)  =  0.0 
sc_by(i,j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  isl,nxl 
do  150  j»l,ny 

if (sc_ax(i, j) .ge.O. )then 
sc_bx (i,j)=do2pl(i,j) 
else 

sc  bx(i , j )=do2pl(i+l, j ) 
endi? 
i50  continue 
c 

do  160  i=l,nx 
do  160  j-l,nyl 
lf(sc_ay(i, j ) .ge.O. ) then 
sc_by(i,j)  .  do2pl(i,j) 


else 

sc_by(i,j)  =  do2pl(i,j+l) 
end  if 
160  continue 
c 
c 

do  170  i»l,nxl 
do  170  j=i,ny 

£o2pl(i,j)  a  sc_ax(i,j)  *  sc  bx(i,i) 
170  continue  ~ 

c 

do  180  i=l,nx 
do  180  j=l,nyl 

go2pl(i,j)  =  sc_ay(i,j)  *  sc  by(i,j) 
180  continue 
c 

do  200  is2,nx 
do  200  j=2,ny 


do2pl(i,j)  a  do2pl(i,j)  -dt*(fo2pl(i,j)-fo2pl(i-l,j))/dx 

'dt*(go2pl(i,j)-go2pl(i,j-l))/dy 
+  +dt*(dopl(i,j)*cl2*dco2(i,j)) 

+  -clt*(do2pl(i,j)*cl6*de(i,j)) 

200  continue 
return 
end 


subroutine  coplusbi 


parameter  (nx*90,ny=20) 
common/gr id/dx , dy , d  t 

common/ fluxC0pi/fcopl(nx,ny) ,gcopl(nx,ny) 
common/dens/dh2  (nx ,  ny  ) ,  dn2  (nx  ,ny  ) ,  dnh3  (nx ,  ny ) , 

$  dh2o(nx,ny),dco(nx,ny),dco2(nx,ny), 

$  dopl ( nx , ny ) , dohpl ( nx , ny ) , dhiopl (nx , ny ) , dnopl ( nx , ny ) , 

S  dnh3pl(nx,ny),dn2pl(nx,ny),de(nx,ny),dh2pl(nx,ny), 

$  do2pl(nx,ny),dcopl(nx,ny),dco2pl(nx,ny),dh3opl(nx,ny), 
S  dphl(nx,ny) i dph2(nx,ny) ,dph3(nx,ny) ,dph4(nx,ny) , 

$  dphS ( nx , ny ) , dph6 ( nx , ny ) , dph7 ( nx , ny ) 
common/veloneu/vxn(nx,ny),vyn(nx,ny) 
common/veloion/vxoxy(nx,ny),vyoxy(nx,ny) 
common/rxra te/cl , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , clO 

common/ rxratebi/cll , cl2 , cl3 , cl4 , cl5 , cl6 , cl7 , cl8 , cl9 
comraon/phrateC0/cph2co,cph3co,cph4co,cph5co 
common/scrat_a/sc_ax(nx,ny),sc_ay(nx,ny) 
common/scrat_b/sc_bx(nx,ny) ,sc  by(nx,ny) 


c -  initialize  scratch  arrays 

c 

do  100  isl^nx 
do  100  j*l,ny 
sc_ax(i,j)  .  0.0 
sc_^ay(i,j)  «  0.0 
sc“bx(i,j)  .  0.0 
sc~by(i,j)  a  0.0 
100  continue 
c 

nxlanx-1 

nyl*ny-l 

c 


130 

140 


150 

c 


160 

c 

c 


170 

c 


180 

c 


& 

& 

& 

& 

& 

200 


c 

c 


$ 

$ 

$ 


do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 
continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 
continue 

do  150  i=l,nxl 
do  150  j=l,ny 
if(sc_axh,j)-ge.0.  )then 
sc_bx(i,j)=dcopl(i, j) 
else 

sc  bx(i,j)=dcopl(i+l,j) 
endil 
continue 

do  160  i=l,nx 
do  160  j»l,nyl 
if(sc_ay(i,j).ge.0.)then 
sc_by(i,j)  =  dcopKi.j) 
else 

sc_by(i,j)  =  dcopl(i,j+l) 
end  if 
continue 


do  170  isl^nxl 
do  170  j*i,ny 

fcopl(i,j)  =  sc  ax(i,j)  *  sc  bx(i,j) 
continue  ■“ 

do  180  i=l,nx 
do  180  j=l,nyl 

gcopl(i,j)  »  sc_ay(i,j)  *  sc  by(i,j) 
continue  “ 


do  200  is2,nx 
do  200  j=2,ny 


dcopl(i,j)  =  dcopl(i,j)  -dt*<£copl(i, j)-fcopl(i- 

-dt*(gcopl(i,j)-gcopl(i,j- 

+dt*(dopl(i,j)*cl7*dco(i,j)) 

+dt*dco(i,j)*((cph2co*dph2(i,j)*(85./110.)) 

+  cph3co*dph3(i, j ) 

+  cph4co*dph4(i,j)  +  cph5co*dph5(i,j)) 
continue 


return 

end 


subroutine  h3oplusbi 

parameter  (nx«90,ny«20) 
common/gr i d / dx , dy , d t 

common/ f luxH30pl/ f h3opl ( nx , ny ) , gh3opl ( nx , ny ) 
common/dens/ dh2 ( nx , ny ) , dn2 ( nx , ny ) , dnh3 ( nx , ny ) , 
dh2o(nx,ny),dco(nx,ny),dco2(nx,ny), 
dopl ( nx , ny ) , dohpl ( nx , ny ) , dh2opl (nx , ny ) , dnopl ( nx , ny 
dnh3pl(nx,ny) ,dn2pl(nx,ny) ,de(nx,ny) ,dh2pl(nx,ny) , 


lj))/dx 

l))/dy 


$  do2pl(nx,ny) ,dcopl(nx,ny) ,dco2pl(nx,ny) ,dh3opl(nx,ny) , 
$  dphl(nx,ny) ,dph2(nx,ny) ,dph3(nx,ny),dph4(nx,ny) , 

$  dph5(nx,ny) ,dph6(nx,ny) ,dph7(nx,ny) 
common/ veloneu/vxn ( nx , ny ) , vyn ( nx , ny ) 
common/ veloi on/ vxoxy ( nx , ny ) , vyoxy ( nx , ny ) 
common/rxrate/cl , c2 , c3 , cA , c5 , c6 , c7 , c8 , c9 , clO 
common/rxratebi/cll,cl2,ci3,cl4,cl5,cl6,cl7 ,cl8,cl9 
common/scrat_a/sc_ax(nx,ny) ,sc_ay(nx,ny) 
common/scrat_b/sc_bx(nx , ny ) , sc_by(nx , ny ) 

c -  initialize  scratch  arrays 

c 

do  100  i=l,nx 
do  100  j=l,ny 
sc_ax(i,j)  =0.0 
sc_ay(i,j)  =0.0 
sc_bx(i,j)  =  0.0 
sc_by(i,j)  =0.0 
100  continue 
c 

nxl=nx-l 

nyl=ny-l 

c 

c 

do  130  i=l,nxl 
do  130  j=l,ny 

sc_ax(i,j)  =  0.5  *  (vxn(i+l,j)  +  vxn(i,j)) 

130  continue 

do  140  i=l,nx 
do  140  j=l,nyl 

sc_ay(i,j)  =  0.5  *  (vyn(i,j+l)  +  vyn(i,j)) 

140  continue 
c 

do  150  i=l,nxl 
do  150  j=l,ny 

if (sc_ax(i, j ) .ge.O. )then 
sc_bx( i , j )=dh3opi ( i , j ) 
else 

sc  bx(i, j)=dh3opl(i+l,j) 
endi? 

150  continue 
c 

do  160  i=l,nx 
do  160  j=l,nyl 
if(sc_ay(i,j).ge.0. )then 
sc_by(i,j)  =  dh3opl(i,j) 
else 

sc_by(i,j)  =  dh3opl(i, j+1) 
endif 
160  continue 
c 
c 

do  170  i=l,nxl 
do  170  j«i,ny 

fh3opl(i,j)  *  sc_ax(i,j)  *  sc_bx(i,j) 

170  continue 

do  180  i=l,nx 
do  180  j.l,nyl 

gh3opl(i,j)  .  sc_ay(i,j)  *  sc_by(i,j) 


c 


180  continue 
c 

do  200  i=2,nx 
do  200  j=2,ny 


dh3opl(i,j)  =  dh3opl(i,j)  -dt*(fh3opl(i,j)-fh3opl(i-l, j))/dx 

-dt*(gh3opl(i,j)-gh3opl(i,j-l))/dy 
+dt*(dh2o(i, j)*cl8*dh2opl(i, j)) 

-dt*(dh3opl(i , j )*cl9*de(i , j ) ) 
continue 
return 
end 


